Adrien Delaye
email: adrien.delaye@irisa.fr

Rafik Sekkal
email: rafik.sekkal@irisa.fr

Eric Anquetil
email: eric.anquetil@irisa.fr

Continuous Marking Menus for Learning Cursive Pen-based Gestures

Keywords: ACM Classification H, 5, 2 Information Interfaces and Presentation: Interaction styles; H, 5, 2 Information Interfaces and Presentation: Input devices and strategies Design, Human Factors

In this paper, we present a new type of Marking menus. Continuous Marking Menus are specifically dedicated to penbased interfaces, and designed to define a set of cursive, realistic handwritten gestures. In menu mode, they offer a continuous visual feedback and fluent exploration of menu hierarchy, inviting the user to execute cursive gestures for invoking the desired commands. In marking mode, a specific gesture recognition method is proposed and proved to be very efficient for recognizing cursive gestures.

INTRODUCTION

Pen-based interfaces offer a very efficient interaction mode, as gestures are realized directly on the target area of the screen, with a stylus manipulated in a natural way. Moreover, gestural commands allow eyes-free selection and are thus specifically suitable for interacting with mobile devices [START_REF] Kurtenbach | The limits of expert performance using hierarchic marking menus[END_REF]. With the progress of techniques for gesture recognition, applications now offer a vocabulary of complex and cursive gestures exploiting all the variety permitted by handwriting, including singularities, loops, inflections points. . . [START_REF] Almaksour | Evolving fuzzy classifiers: Application to incremental learning of handwritten gesture recognition systems[END_REF][START_REF] Rubine | Specifying gestures by example[END_REF][START_REF] Willems | Iconic and multi-stroke gesture recognition[END_REF]. Such richer sets of gestures allow higher symbolic significance, and can be easier to remember than piecewise linear gestures. Despite the aforementioned advantages, gesturebased interfaces suffer that gestures are not self-revealing (list of possible gestures is not easily accessible to the user). To cope with this, shapes of gestures are often chosen to be as much as possible linked by a logical relation with the command to be executed (semantic, visual or mnemonic relation). For instance, the Copy action can be associated to c ACM, 2011. This is the author's version of the work. It is posted here by permission of ACM for your personal use. Not for redistribution. The definitive version was published in Proceedings of the 15th International Conference on Intelligent User Interfaces, http://doi.acm.org /10.1145/1943403.1943455 a "C"-shaped gesture, or Erase action to a "cross-out" gesture. The more an action-gesture association respects some relation understandable for the user, the easier it is to learn and remember it [START_REF] Rubine | Specifying gestures by example[END_REF][START_REF] Willems | Iconic and multi-stroke gesture recognition[END_REF]. However, this logical association is often insufficient, and learning or remembering gestures can remain a tedious task, especially if the number of gestures is important. The logic of action-gestures associations can appear as poorly intuitive for some users, and is often unable to deal with purely abstract or non-representable actions. Alternative solutions exist for making gestures self-revealing. In [START_REF] Bau | OctoPocus: a dynamic guide for learning gesture-based command sets[END_REF], authors propose an interactive guiding of pen trajectory by dynamically displaying gestures suggestions. This assistance facilitates memorization and no effort is needed to remember gestures. Another solution is to let users define their own gestures and associate actions, but this requires recognition engines able to deal with new shapes, while demanding as few training samples as possible. Designing this type of recognition methods is still a challenging pattern recognition problem [START_REF] Almaksour | Evolving fuzzy classifiers: Application to incremental learning of handwritten gesture recognition systems[END_REF].

Marking menus offer an alternative for assisting the user in learning gestures [START_REF] Kurtenbach | The limits of expert performance using hierarchic marking menus[END_REF]. These menus offer two modes : item selection from a menu, and gesture command without the menu being displayed. First mode is a traditional use of a pie menu: user selects an item with a pen (or a sub-menu he wants to explore in the case of hierarchical menus). Second mode implies gesture recognition, and works like a gesturebased interface, where the user is familiar with the gestures commands for triggering actions. Kurtenbach's idea is to link the gesture shapes to the layout of the underlying menu: the physical movement for selecting a menu item is identical to the physical movement required to make the associated gesture command [START_REF] Kurtenbach | The limits of expert performance using hierarchic marking menus[END_REF]. By repetitively using the menu mode, user can thus memorize gestures without explicitly making any effort. One drawback of marking menus and variants, though, is they are limited in expressiveness as they are restricted to polygonal gestures [START_REF] Kurtenbach | The limits of expert performance using hierarchic marking menus[END_REF], or sequences of linear gestures [START_REF] Zhao | Simple vs. compound mark hierarchical marking menus[END_REF][START_REF] Bailly | Wave menus: improving the novice mode of hierarchical marking menus[END_REF] (except for Flower Menus that allow degrees of curvature [START_REF] Bailly | Flower menus: a new type of marking menu with large menu breadth, within groups and efficient expert mode memorization[END_REF]).

In this paper, we introduce Continuous Marking Menus, an improved type of marking menus implying a rich set of cursive gestures, that can be executed in a single continuous gesture with a pen or finger. We first propose some novelties for topology and interaction behavior to help the user make selection of menu items in a more fluent and gesturelike fashion. Integration of specific inertia areas permits to absorb imprecision in the change of direction, coupled with a continuous visual feedback, supports the user in the making of fluent gestures. Finally, we propose a specific gesture recognition technique that can handle the variability of resulting curved gestures in the marking mode. A review of marking menus related work is given in the first section, before we introduce Continuous Marking Menus design. The recognition method used in the marking mode is then presented briefly, together with an experiment that show its ability to properly recognize cursive gestures.

RELATED WORK

Marking Menus were introduced for making gestures selfrevealing with the help of a radial menu that is popped-out according to the user need. They permit to implicitly learn the gesture commands by repetitively using the menu mode, and the transition to marking mode is facilitated by making user rehearse the physical movement involved in making the gesture every time a selection from the menu is made [START_REF] Kurtenbach | The limits of expert performance using hierarchic marking menus[END_REF]. The basic, non-hierarchical version of marking menus involves simple radial gestures from menu center toward the chosen direction. In the hierarchical version (Compound Marking Menus), several levels can be explored by making a piecewise linear gesture, each part corresponding to the selection of an item or sub-menu at one level. Such gestures can be easily executed with the help of a mouse, even if the use of the pen is significantly more efficient [START_REF] Kurtenbach | The limits of expert performance using hierarchic marking menus[END_REF]. However, strong constraints about their shapes (polygonal shape with abrupt angles) are actually in contradiction with the drawing of fluent handwritten gestures with a pen, because users have to pause at each menu level, and gestures are more a sequence of linear movements. Simple Marking Menus [START_REF] Zhao | Simple vs. compound mark hierarchical marking menus[END_REF], or variants as Zone and Polygon menus [START_REF] Zhao | Zone and polygon menus: using relative position to increase the breadth of multi-stroke marking menus[END_REF] offer to decompose the selection of sub-menus and items into several disconnected marks, one for each level. They permit to solve two problems of Compound Marking Menus. First, they reduce the space required for menu display and drawing of gestures in both modes. Secondly, they allow an increased depth (more successive levels in the hierarchy), by simply adjusting the number of menu selections (or marks) in the sequence. Zone and Polygon menus further permit more items at each level (increase breadth). Simple Marking Menus and alike, however, are less efficient in menu mode because the navigation path through the hierarchy is not visible to the user, and no previsualization (browsing of sub-menus for a fast inspection of available commands) is possible. More importantly, by nature they require several distinct gestures to invoke commands, which is a step away from our objective of permitting fully continuous gestures. Wave Menus [START_REF] Bailly | Wave menus: improving the novice mode of hierarchical marking menus[END_REF] As a summary, the existing variants of Marking Menus only offer a rather limited vocabulary of gestures, mostly linear or piecewise linear gestures, not dedicated to direct interaction with pen or finger. Our work aims at improving the marking mode of these menus, by tolerating a much larger range of gestures expected to be faster and more efficient to execute. The challenge is twofold : how to design a new type of marking menu offering an enlarged gestural expressiveness while maintaining an efficient menu mode, and how this menu mode can be optimized for inviting to make fluent, gestural-like menu selections.

CONTINUOUS MARKING MENUS

Continuous Marking Menus are a variant of Compound Marking Menus [START_REF] Kurtenbach | The limits of expert performance using hierarchic marking menus[END_REF] with enhanced interaction and continuous user feedback in the menu mode that facilitates continuous and fluent movements for selecting commands.

Design

Continuous Marking Menus are based on hierarchical radial menus with a breadth of 4, 6 or 8 items at each level. All the examples shown in this article are based on a 6-breadth version, which we used for experiments. For each level, the menu consists of a set of kite-shaped branches, centered around the point where the menu was triggered by the cursor. The labels are displayed at the extremity of each branch, in a rectangle oriented perpendicularly to the direction of the branch. Contrarily to traditional pie menus, the active areas are not contiguous (there is an inactive blank zone between two neighboring branches), in order to guide the user make menu selections in a more trajectory-like way. While classical Marking Menus display only labels as targets to be selected, we propose to display directed branches that suggest a path to invoke a command. Figure 1 shows two configurations of this menu, with 6 branches each. The two configurations only differ by a rotation (axis-shift), and they will be alternately displayed when exploring the menu hierarchy levels, in order to avoid ambiguity problems by forcing a change of direction whenever a sub-menu is developed.

Interaction

A specific interactive behavior is required for assisting the user and invite him to make continuous gestures while exploring the menu. Specifically, Continuous Menus offer a Pre-selection area is useful for browsing the menu level and let the user look for desired item or sub-menu. The associated visual feedback highlights the branch more and more intensely as the user cursor moves toward the branch extremity. Simultaneously, the other branches of the menu level are faded out, as shown in figure 2(b). When the cursor crosses pre-selection limit, branch menu is considered as pre-selected, and the other alternative branches are completely invisible. As the cursor enters the inertia area, the sub-menu of the pre-selected branch fades in and its branches are more and more visible as the cursor moves away from the center of the inertia area (see figure 2(c)). In this area, the center of the sub-menu sticks to the movement of the cursor, until it crosses the outter boundaries of the area (see figure 2(d)). This area permits to absorb the inertia of the pen movement when a change of direction is executed in order to encourage the user to make fluent moves when exploring the menu hierarchy. Indeed, thanks to this mechanism, no stopping or dwelling in the user movement is needed when changing direction. Inertia area can be seen as the tolerance area where the change of direction can be done.

Detailed menu selection process

The dynamic interaction offered by the menu when user make a menu selection can be described in details by the following succession of steps, illustrated in figure 3 1 .

(a) Initial menu display, centered on the cursor. (b) Pre-selection area. The chosen branch ("Europe") is highlighted, the other ones fade out.

1 See video at http://www.irisa.fr/imadoc/CMM.html (c) Pre-selection of sub-menu ("Europe"), and entering in the inertia area. Only the chosen branch is visible. (d) Inertia area. Sub-menu branches progressively appear, and the sub-menu follows the position of the cursor. (e) Crossing of the inertia area limit. The sub-menu is now fixed. The parent branch ("Europe") is validated, and its shape changes to a smaller kite as a visual confirmation.

Another cycle starts (back to step (a)).

At any time before Branch validation (step (e)), the user can move back to the parent menu and change the preselected branch by moving back the cursor along the branches toward the center of parent menu. The continuous feedback, including fading out of parent overlooked branches and fading in of subsequent child sub-menu not only stimulates the user for pursuing his gesture, but is also precious for maintaining an acceptable visual complexity. The chart in figure 3 shows that the number of items displayed during the selection process remains limited, and branches from two level of menus are never displayed simultaneously.

MARKING MODE WITH CURVED GESTURES

The marking mode of Continuous menus involves a gesture recognizer that do not provide any visual interaction during the execution of the gesture. Once the gesture is accomplished, several pre-processing are applied (normalization, resampling), and the Dynamic Type Warping (DTW) recognition technique is used : elastic matching distance is computed by dynamic programming between the input gesture and each gesture prototype (modeling the ideal gesture for each possible command). The input gesture is assigned the class corresponding to the closest prototype according to DTW distance, and the associated command is invoked. Since elastic matching distance does not involve a polygonal approximation of gestures, it is suited to handle properly cursive and distorted gestures. Moreover, the nearest-neighbor classification paradigm require no training phase, so it could be possible for the user to dynamically modify, add or remove some menu items according to his needs. A database of gestures was constituted by asking subjects to execute 32 different gestures associated to a 3-level menu with a breadth of 6 items per level (see figure 4(b)). In order to assist participants in achieving expert behavior, we displayed a visual guide for the required mark, thus eliminating the need for familiarity with gestures. The visual guide corresponds to the display of starting and ending point and several fuzzy grey areas denoting the inertia areas, where the user should approximately make direction changes (figure 4(a) shows an example of visual guide). A Tablet PC with a pen was used for the input. 10 people participated in the experiment, with age between 18 and 45. Each of them executed 10 instances of each gesture, amounting to a total number of 3200 samples. Our experiment with DTW recognizer reached a correct recognition rate of 99,97%, confirming its great ability to discriminate between cursive gestures.

CONCLUSION

In this work we introduced Continuous Marking Menus for improving the expressiveness of marking menus and admitting a vocabulary of cursive, fluent handwritten gestures, and thus offer a more natural and efficient interaction in the marking mode. We proposed an improved interaction in the menu mode and a specific mechanism for handling imprecise direction changes, suited to the specificities of pen based or finger based gestures. The method proposed for recognition of gestures in the marking mode is well adapted to cursive gestures, and its performance was experimentally validated. The very high recognition rate suggests that the breadth and depth choices for our menus are not a limitation for our classification method. Future work will focus on evaluating acceptability and efficiency of these menus in comparison with alternative techniques. The gain of efficiency and user satisfaction should be measured and compared in the marking mode when using cursive vs polygonal gestures. Finally, in order to push further the idea of efficient single-gesture commands, we will also consider combining both the selection and control of actions in a single gesture, in the light of the works of Pook [START_REF] Pook | Control menus: excecution and control in a single interactor[END_REF].

Figure 1 .

 1 Figure 1. Two configurations of Continuous Menus alternately used in successive hierarchical levels: level i (a) and sub-menu level i+1 (b)

Figure 2 .

 2 Figure 2. Interaction areas

Figure 3 .

 3 Figure 3. Steps for menu selection (a)-(f) and variation of the visual complexity.

Figure 4 .

 4 Figure 4. Example of guide for assisting expert behavior (a), and table of 32 gestures classes (b).

ACKNOWLEDGMENTS

This work was partly supported by the Region Bretagne through the ScriptEverywhere project.


Brunella Pilato

Simona De Summa

Katia Danza

Rossana Lambo

Angelo Paradiso

PhD Stefania Tommasi
email: s.tommasi@oncologico.bari.it

MATERNAL AND PATERNAL LINEAGE DOUBLE HETEROZYGOSITY ALTERATION IN FAMILIAL BREAST CANCER: A FIRST CASE REPORT

Keywords: BRCA1, BRCA2, Mutation, Breast and Ovarian cancer, Double Heterozygosity (DH)

Hereditary breast cancer syndrome was firstly associated with BRCA1 and BRCA2 genes whose mutations confer high risk to develop breast and/or ovarian cancer. Double heterozygosity is a rare condition in which both BRCA1 and BRCA2 mutations are present in a family at the same time.

In the current study, a family with double heterozygosity has been reported. Furthermore, for the first time a molecular analysis in both proband lineages, maternal and paternal, has been reported to understand the provenience of both germinal mutations.The case regards a woman that developed breast and ovarian cancer with liver metastasis which presented two mutations, each in the two genes, transmitted from her mother and her father respectively. In this family all available members have been investigated. The concomitantly presence of these peculiar mutations was never reported before suggesting a link with Caucasian population from Southern Italy.

Introduction

High penetrating BRCA1 and BRCA2 genes germline mutations have been associated with hereditary breast cancer risk. Women with germline BRCA mutations have a relevant risk increase for breast and ovarian cancer if compared with the global population. It has been estimated that about 10-15% of cancer patients have mutations in BRCA1 and 7% in BRCA2 genes [START_REF] Tommasi | Molecular and in silico analysis of brca1 and brca2 variants[END_REF]. There are different types of mutations listed in Breast Cancer Information Core professional databases (BIC).

Most are small insertions or deletions generating truncated proteins with missing common functionalities. Moreover, missense mutations, intronic variants sequences (IVS) and unclassified variants (UV) have been identified, reporting a not clear role in breast cancer susceptibility [START_REF] Diez | Identification of a new complex deleterious mutation in exon 18 of the brca2 gene in a hereditary male/female breast cancer family[END_REF].

Previous studies has highlighted that 80% of patients with a cancer family history and a high BRCA mutation risk have at least one polymorphism that could increase breast cancer risk [START_REF] Tommasi | Molecular and in silico analysis of brca1 and brca2 variants[END_REF][START_REF] Pilato | Mutations and polymorphic brca variants transmission in breast cancer familial members[END_REF]. Furthermore, screening the mutations in the two genes and approaching with new technical methods, large genomic rearrangements have been identified. In particular, at least 81 different large genomic rearrangements (LGRs) in BRCA1 gene and few rearrangements in BRCA2 gene have been characterized [START_REF] Sluiter | Large genomic rearrangements of the brca1 and brca2 genes: Review of the literature and report of a novel brca1 mutation[END_REF].

However, the probability of having both genes mutated in the same individual is very low and only few cases showing double heterozigosity (DH) in Italian population have been described [START_REF] Zuradelli | Four new cases of double heterozygosity for brca1 and brca2 gene mutations: Clinical, pathological, and family characteristics[END_REF][START_REF] Cortesi | Italian family with two independent mutations: 3358t/a in brca1 and 8756dela in brca2 genes[END_REF][START_REF] Musolino | A breast cancer patient from italy with germline mutations in both the brca1 and brca2 genes[END_REF].

In the current study, one case of breast cancer's affected family with an individual with DH for BRCA1 and BRCA2 gene mutations has been reported.

Methods

The proband and related family members have been enrolled from National Cancer Centre of Bari -Counseling Program and all of them undersigned an informed consent, previously approved by the ethical committee, to perform molecular analyses and to use data for research purposes [START_REF] Tommasi | Brca1 mutations and polymorphisms in a hospital-based consecutive series of breast cancer patients from apulia, italy[END_REF]. The patient was also characterized according to pathological features [START_REF] Berry | Brcapro validation, sensitivity of genetic testing of brca1/brca2, and prevalence of other breast cancer susceptibility genes[END_REF]; this is a smoker patient (5 cig/die). DNA was extracted using QIAmp DNA blood midi kit (Qiagen, Valencia, CA, USA) following manufacturer's instructions. Molecular analysis were performed as previously described [START_REF] Tommasi | Molecular and in silico analysis of brca1 and brca2 variants[END_REF][START_REF] Diez | Identification of a new complex deleterious mutation in exon 18 of the brca2 gene in a hereditary male/female breast cancer family[END_REF][START_REF] Tommasi | Brca1 mutations and polymorphisms in a hospital-based consecutive series of breast cancer patients from apulia, italy[END_REF]].

Family's information has allowed to design a pedigree with the Progeny v6 Program.

The variants in the sequence were characterized and compared to those present in the online databases BIC [http://www.research.nhgri.nih.gov/bic/].

Case Report

This case is a 44 years old female proband that was enrolled in the Genetic Counseling Program of NCC Bari during 2008. In 2004, the patient has undergone to a nodulectomy on the left side and in 2005 to a quadrantectomy. Histopathological analysis showed a medium grade intraductal carcinoma with negative ER, PgR and HER-2. Subsequently, she underwent a local radiotherapy. In 2008, she was affected by a third grade bilateral papillary ovarian adenocarcinoma. Analyzing peritoneal washing liquid, ER and PgR values were negative. In 2010, metastases of ovarian origin have been diagnosed in the liver. Molecular analysis showed a double heterozygosity for BRCA1 and BRCA2 genes: the BRCA1 NT_5382insC and the BRCA2 NT_6024delTA mutations. During the counseling interview, the patient referred a maternal history of breast and ovarian cancer. Her mother developed breast and ovarian cancer at 58 and 60 years respectively. For this family, the proband's mother, father, one unaffected sister and two unaffected brothers were available for genetic analysis. Both mutations resulted to be inherited, one from the mother (BRCA1 NT_5382insC) and one from the father (BRCA2 NT_6024delTA). With exception of the proband, that had a DH for BRCA1 and BRCA2 gene mutations, only the unaffected proband's sister presented BRCA2 NT_6024delTA mutation. Figure 1 shows the genealogic tree of this family.

Discussion

In the current study, one case of DH (BRCA1 NT_5382insC and BRCA2 NT_6024delTA) in a patient belonging to family with breast and ovarian cancer history has been reported. Interestingly, analyzing both maternal and paternal lineages we evidenced that the two germinal mutations were respectively inherited one by the mother and one by the father. The double heterozygosity is the presence of pathological mutations in both BRCA1 and BRCA2 genes in the same subject and it is an extremely rare event mostly shown in Ashkenazi Jewish [START_REF] Leegte | Phenotypic expression of double heterozygosity for brca1 and brca2 germline mutations[END_REF]. The BRCA1 NT_5382insC is the most common mutation in Caucasian population from Puglia [START_REF] Tommasi | Molecular and in silico analysis of brca1 and brca2 variants[END_REF] and in Askenazi in which accounted for 26 (10.1%) of all the BRCA1 founder mutations, appearing in 60% of ovarian cancer and 30% of breast cancer in Ashkenazi population [START_REF] Abeliovich | The founder mutations 185delag and 5382insc in brca1 and 6174delt in brca2 appear in 60% of ovarian cancer and 30% of early-onset breast cancer patients among ashkenazi women[END_REF]. The BRCA1 NT_5382insC mutation is located near the 3' end of the BRCA1 coding region and consists of an insertion of the nucleotide deoxycytidylate (C) at nucleotide position 5382 of the gene. Analyzing a possible correlation between this BRCA1 mutation and specific phenotypes, Rennert et al. [START_REF] Rennert | Differences in the characteristics of families with brca1 and brca2 mutations in israel[END_REF] showed as the BRCA1 NT_5382insC mutation was more expressed in bilateral breast cancer (38.9% of families, 15.4% of women with cancer, 6.1% of all women in family) and metachronous breast and ovary tumors (22.2%, 9.8% and 4.5% correspondingly). Moreover, studies about ovarian cancer's lifetime risk and inherited BRCA1 NT_5382insC mutation in family members and in population series reported a lifetime risk of ovarian cancer of 28-66% [START_REF] Satagopan | Ovarian cancer risk in ashkenazi jewish carriers of brca1 and brca2 mutations[END_REF] The second alteration found in this studied proband was BRCA2 NT_6024delTA mutation (Accession Number BIC database: 4752). This mutation consists of a deletion of the nucleotides TA that generates the truncation of the BRCA2 protein at codon 1943 [START_REF] Kataki | Identification of germline brca1 and brca2 genetic alterations in greek breast cancer moderate-risk and low-risk individuals--correlation with clinicopathological data[END_REF]. It was firstly identified in Italian Breast and/or Ovarian Carcinoma patients in 1998 [START_REF] Benedetti | Characterization of ten novel and 13 recurring brca1 and brca2 germline mutations in italian breast and/or ovarian carcinoma patients. Mutations in brief no. 178[END_REF], but subsequently it has been detected in the ovarian cancer cluster region (OCCR) of Greek patients [START_REF] Armakolas | Brca2 gene mutations in greek patients with familial breast cancer[END_REF], and in other different population with breast cancer [START_REF] Tommasi | Molecular and in silico analysis of brca1 and brca2 variants[END_REF]. No published data demonstrated a founder effect of this BRCA2 mutation.

The current study is the first describing these two mutations in DH and it is the first registration of a DH case in Southern Italian population.

Even if rare, most common mutations reported in DH were BRCA1 NT_185delAG and BRCA2 NT_6174delT [START_REF] Moslehi | An unaffected individual from a breast/ovarian cancer family with germline mutations in both brca1 and brca2[END_REF] which presented a typical phenotype [START_REF] Leegte | Phenotypic expression of double heterozygosity for brca1 and brca2 germline mutations[END_REF][START_REF] Moslehi | An unaffected individual from a breast/ovarian cancer family with germline mutations in both brca1 and brca2[END_REF]. Indeed, BRCA1 NT_5382insC has been reported in DH very few times [START_REF] Leegte | Phenotypic expression of double heterozygosity for brca1 and brca2 germline mutations[END_REF][START_REF] Bell | Selective loss of heterozygosity in multiple breast cancers from a carrier of mutations in both brca1 and brca2[END_REF] showing biallelic inactivation of either BRCA1 or BRCA2 thus suggesting that these genetic events are functionally equivalent in initiating tumorigenesis.

By reviewing published data about Ashkenazi Jewish descent, Leegte et al. [START_REF] Leegte | Phenotypic expression of double heterozygosity for brca1 and brca2 germline mutations[END_REF] analyzes 4 cases of double heterozygosity and reviews 30 cases reported in literature in order to assess phenotypic expression with respect to age of onset, cumulative lifetime risk and multiple primary tumors. The phenotypic expression in this small series of 34 female cases of DH appeared to be highly variable, however, analyzing the cancer incidence in DH at age 70 he reported a value of 80% for breast cancer and 84% for breast and ovarian. Struewing et al [19] reported a risk of 40-73% at age 70 for BRCA1. Warner et al [START_REF] Warner | Prevalence and penetrance of brca1 and brca2 gene mutations in unselected ashkenazi jewish women with breast cancer[END_REF] calculate a risk of 60% at age 70 for BRCA1 while Satagopan et al 31-80% at age 70 for BRCA1 and 14-50% for BRCA2 [START_REF] Satagopan | The lifetime risks of breast cancer in ashkenazi jewish carriers of brca1 and brca2 mutations[END_REF]. These data were comparable to the family based data of Brose et al [START_REF] Brose | Cancer risk estimates for brca1 mutation carriers identified in a risk evaluation program[END_REF], with 37% Ashkenazi cases (74-82% by age 70 for BRCA1), and the population based data of King et al (69% for BRCA1 and 74% for BRCA2 at age 70) [START_REF] King | Breast and ovarian cancer risks due to inherited mutations in brca1 and brca2[END_REF].

Only few published data concerned DH incidence in non-Ashkenazi and in Italian population. Claus et al [START_REF] Claus | Prevalence of brca1 and brca2 mutations in women diagnosed with ductal carcinoma in situ[END_REF] reported one case of double heterozygosity in both BRCA1 and BRCA2 genes (BRCA1 W321X and BRCA2 NT_3398del5 mutations) and a woman that had 2 mutations in BRCA2: BRCA2 S1630X and BRCA2 Q1782K. Regarding the double heterozygosity in Italian population, Cortesi et al [START_REF] Cortesi | Italian family with two independent mutations: 3358t/a in brca1 and 8756dela in brca2 genes[END_REF] studied an Italian family where two different non-Ashkenazi Jewish mutations were present: 3358T/A in BRCA1 and NT_8756delA in BRCA2. The BRCA1 mutation was inherited since 50% of the members studied were carriers whereas the BRCA2 NT_8756delA mutation seemed to be a "de novo" mutation. Musolino et al [START_REF] Musolino | A breast cancer patient from italy with germline mutations in both the brca1 and brca2 genes[END_REF] studied a breast cancer patient from Italy with BRCA1 Y1429X and BRCA2 Q2580X mutations reporting the first case in Italy of a non-Ashkenazi double heterozygosity for BRCA1 and BRCA2 genes. At last, Zurardelli et al [START_REF] Zuradelli | Four new cases of double heterozygosity for brca1 and brca2 gene mutations: Clinical, pathological, and family characteristics[END_REF] discussed four new cases of double heterozygosity showing a non correlation of age of onset, specific phenotype and a presence of pathological mutations in both genes at the same time. The studied patient also presented a severe phenotype developing breast and ovarian cancer at 38 and 42 years respectively, both characterized by having absence of hormone receptor and low tumor differentiation grade. Furthermore, she had a bilateral ovarian carcinoma and discovered to have liver metastases of ovarian origin at 44 years. Moreover, the concomitant presence of BRCA1 NT_5382insC and the BRCA2 NT_6024delTA could predispose to a high incidence of ovarian cancer.

In conclusion, our report confirmed one more time the BRCA mutation type dependency on ethnicity. Furthermore, studying BRCA mutations incidence in breast cancer families, in particular when gene mutations are present in DH has both scientific and clinical relevance. From the scientific point of view it is important to understand the interaction between BRCA1 and BRCA2 which participate the same metabolic pathway in damage repair. From the clinical point of view, DH patients should undergo a peculiar follow-up and surveillance strategies and, in particular, because most of the DH cases presented in the literature were affected by both breast and ovary cancers, a prophylactic oophorectomy could be suggested.

Figure legend Figure 1 .

 legend1 Figure legend

Ackowledgements

Partially supported by Regione Puglia (DIEF 2007)


Vandna Jowaheer

Brajendra Sutradhar
email: bsutradh@math.mun.ca

GMM versus GQL inferences for panel count data

Keywords: Efficiency, Estimation methods, Gaussian random effects, Nonstationary longitudinal correlation structure, Repeated count responses

It is well known that the likelihood inferences in dynamic mixed models for count data is extremely complicated. In this paper, we, first, develop a generalized method of moments (GMM) approach for the estimation of the parameters of such models. We then consider an alternative generalized quasi-likelihood (GQL) approach.

The relative efficiency of the GQL approach to the GMM approach is examined by comparing the asymptotic variances of the GQL estimates of the parameters to the corresponding asymptotic variances of the GMM estimates.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

INTRODUCTION

In longitudinal studies for count data, a small number of repeated count responses along with a set of covariates are collected from a large number of individuals. For example, in health care utilization data, the number of visits to the physician by a large number of independent individuals may be recorded over a period of several years. Also, the information on the covariates-gender, education level and age may be recorded for each individual. In this set up, it is most likely that the repeated count responses of an individual may be influenced by individual's random effect and conditional on such a random effect the repeated responses are auto-correlated. This type of panel or longitudinal count data have been studied by many authors such as Sutradhar and [START_REF] Sutradhar | On Generalized Quasilikelihood Inference in Longitudinal Mixed Model for Count Data[END_REF] [see also [START_REF] Wooldridge | Distribution-free estimation of some non-linear panel data models[END_REF]]. In notation, let (y i1 , . . . , y it , . . . , y iT ) denote the T repeated count responses for the ith subject, i = 1, . . . , K. Also, let x it be the p × 1 vector of covariates corresponding to y it , and β is the p × 1 regression effects of x it on y it . Next suppose that in addition to x it , the repeated responses of the ith individual are also influenced by one random effect γ * i . Assume that γ * i iid ∼ N (0, σ 2 γ ). Now, by writing the marginal density of y i1 (but conditional on γ * ) as f i1 (y i1 |γ * i ), and the conditional density of y it given y i,t-1 as f it|t-1 (y it |y i,t-1 , γ * i ) for t = 2, . . . , T, one may write the likelihood function for the ith individual as

L i (β, ρ, σ γ ) = ∞ -∞ f i1 (y i1 |γ * i )Π T t=2 f it|t-1 (y it |y i,t-1 , γ * i ) × φ(σ γ γ i )dγ i , (1.1)
φ(γ i ) being the standard normal density, and ρ being the dependence parameter of y it on y i,t-1 conditional on γ * i , yielding the likelihood function

L(β, ρ, σ γ |y) = Π K i=1 L i (β, ρ, σ γ ) (1.2)
for the complete data set y.

Note that since y i1 conditional on γ i = γ * i /σ γ has the Poisson density, we write

f i1 (y i1 |γ i ) = exp(-µ * i1 )µ * i1 y i1 y i1 ! , (1.3) where µ * i1 = exp(x i1 β + σ γ γ i ).
Next suppose that for t = 2, . . . , T, y it and y i,t-1 has the dynamic relationship

[y it |γ i ] = ρ • [y i,t-1 |γ i ] + [d it |γ i ],
(1.4) [START_REF] Mckenzie | Some ARMA Models for Dependent Sequences of Poisson Counts[END_REF], [START_REF] Sutradhar | An Review on Regression Models for Discrete Longitudinal Responses[END_REF], [START_REF] Sutradhar | On Generalized Quasilikelihood Inference in Longitudinal Mixed Model for Count Data[END_REF]] where ρ

• y i,t-1 = y i,t-1 s=1 b s (ρ) with P r[b s (ρ) = 1] = ρ and P r[b s (ρ) = 0] = 1 -ρ. In (1.3) [y i1 |γ i ] ∼ P oi(µ * i1 ), and in (1.4) for t = 2, . . . , T, [d it |γ i ] ∼ P oi(µ * it -ρµ * i,t-1 ), with µ * it = exp(x it β + σ γ γ i ), for t = 1, . . . , T.
(1.5)

By (1.4), it can be shown that conditional on γ i for t = 2, . . . , T, the lag 1 conditional density has the form

f it|t-1 (y it |y i,t-1 , γ i ) = exp[-(µ * it -ρµ * i,t-1 )] × min(y it ,y i,t-1 ) s it =0 y i,t-1 !ρ s it (1 -ρ) y it -s it (µ * it -ρµ * i,t-1 ) y it -s it s it !(y i,t-1 -s it )!(y it -s it )!
, (1.6) [START_REF] Freeland | Analysis of Low Count Time Series Data by Poisson Autoregression[END_REF] where µ * it = exp(x it β + σ γ γ i ) for all i = 1, . . . , K; t = 1, . . . , T. It is clear that when this conditional density (1.6) is used in (1.1), it provides a complicated likelihood function by (1.2). Thus it is impractical to use the likelihood approach for such repeated count data analysis.

Note that even though there may be other ways to model the repeated count data, the auto-regressive order 1 (AR(1)) type dynamic model (1.4) is much more attractive in practice because in stationary case it yields a correlation pattern where the lag correlations decay exponentially as lag increases. This is one of the main reasons why this model is used frequently [START_REF] Mckenzie | Some ARMA Models for Dependent Sequences of Poisson Counts[END_REF], [START_REF] Freeland | Analysis of Low Count Time Series Data by Poisson Autoregression[END_REF], [START_REF] Sutradhar | On Generalized Quasilikelihood Inference in Longitudinal Mixed Model for Count Data[END_REF]) as compared to other such as moving average order 1 (MA(1)) or equi-correlations type models. But, because of the complexity of the likelihood approach, it is of interest to seek for alternative approaches for the estimation of the parameters of the dynamic model (1.4).

With regard to the estimation of the parameters in a class of linear correlated dynamic mixed models, many econometricians such as [START_REF] Arellano | Some Tests of Specification for Panel Data : Monte Carlo Evidence and an Application to Employment Equations[END_REF], [START_REF] Ahn | Efficient Estimation of Models for Dynamic Panel Data[END_REF], [START_REF] Blundell | Initial Conditions and Moment Restrictions in Dynamic Panel Data Models[END_REF][START_REF] Blundell | Initial Conditions and Moment Restrictions in Dynamic Panel Data Models[END_REF][START_REF] Imbens | Generalized Method of Moments and Empirical Likelihood[END_REF] have exploited the well known generalized method of moments (GMM) due to [START_REF] Hansen | Large-Sample Properties of Generalized Method of Moment Estimators[END_REF]. But there does not appear any discussion about the use of this GMM approach in the present non-linear correlated dynamic mixed model set up for the count data. In this paper, we, first, develop this GMM approach to deal with estimation of the parameters of the proposed non-linear model for the repeated count data. Next, we consider a recently developed generalized quasi-likelihood (GQL) approach [START_REF] Sutradhar | On Generalized Quasilikelihood Inference in Longitudinal Mixed Model for Count Data[END_REF] and [START_REF] Sutradhar | An Review on Regression Models for Discrete Longitudinal Responses[END_REF]) for the inferences for correlated data. Details on the GMM and GQL estimation approaches for panel count data are given in Section 2. Note that the main objective of the paper is to examine the relative efficiency of the GQL and GMM approaches for the panel count data model (1.4). We also note here that recently [START_REF] Sutradhar | Generalized method of moments versus generalized inferences in binary panel data models[END_REF] have examined the relative efficiency of the GMM and GQL approaches in a binary panel data set up. But, as it is difficult to make a mathematical connection between the panel count model (1.4) and the binary panel data model considered by [START_REF] Sutradhar | Generalized method of moments versus generalized inferences in binary panel data models[END_REF], nothing can be said on efficiency performances under the count data model based on the results for a binary panel data model. Thus, it is important to conduct a separate study for the panel count data models. The efficiency comparisons between the GMM and GQL approaches for the panel count data model (1.4) are given in Section 3.

GMM AND GQL ESTIMATION

Estimation of β

To estimate β, we assume that the other two parameters, namely σ 2 γ and ρ are known. Note that both the GMM and GQL are moments based approaches. Estimation of the β vector requires the first and second order moments of the data under the model.

We first provide these moments and then use them to form the GMM and GQL estimating equations.

First and second order unconditional moments

Let µ it and σ itt denote the mean and the variance of y it . Since γ i i.i.d.

∼ N (0, 1), and because it follows from (1.4) that y it conditional on γ i follows the Poisson distribution with parameter µ * it given in (1.5), y it then unconditionally has the mean and the variance given by,

µ it = E[Y it ] = E γ i E[Y it |γ i ] = E γ i [µ * it ]) = exp[x it β + σ 2 γ /2] (2.1) σ itt = var[Y it ] = E[Y 2 it ] -[E(Y it )] 2 = E γ i E[Y 2 it |γ i ] -µ 2 it = E γ i [µ * it + µ * it 2 ] -µ 2 it = µ it + [exp(σ 2 γ ) -1]µ 2 it (2.2) A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Furthermore, it may be shown from (1.4) that conditional on γ i , the covariance between y iu and y it for (u < t) is given by

cov(Y iu , Y it |γ i ) = ρ t-u µ * iu . (2.3)
Consequently, under the serial dependence model, the unconditional covariance between y iu and y it has the form

σ iut (β, σ 2 γ , ρ) = E[cov{(Y iu , Y it )|γ i }] + cov[E(Y iu |γ i ), E(Y it |γ i )] = ρ t-u µ iu + [exp(σ 2 γ ) -1]µ iu µ it .
(2.4)

GMM estimation for β

Note that when the traditional method of moments (MM) is used to estimate the β vector, one solves the unbiased moment estimating equation

K i=1 ψ i1 (β, σ 2 γ ) = 0, (2.5)
where [START_REF] Jiang | Consistent Estimators in Generalized Linear Mixed Models[END_REF], [START_REF] Sutradhar | On Exact Quasilikelihood Inference in Generalized Linear Mixed Models[END_REF]] is an unbiased moment function as

ψ 1i (β, σ 2 γ ) = T t=1 [x it (y it -µ it )] [
E[Y it ] = µ it leading to E[ψ 1i (β, σ 2 γ )] = 0.
In notation of (2.5) but for known σ 2 γ , in the GMM approach [START_REF] Hansen | Large-Sample Properties of Generalized Method of Moment Estimators[END_REF]) one would minimize the distance function

Q(β) = K -1 K i=1 ψ i1 (β|σ 2 γ )) C 1 K i=1 ψ i1 (β|σ 2 γ ) , (2.6)
where C 1 is a suitable weight. An optimal choice for C 1 would be the inverse of the variance of the unbiased moment function, that is,

C 1 = [var(K -1 K i=1 ψ i1 (β, σ 2 γ ))] -1 . (2.7)
Note that in the context of similar longitudinal studies, some authors such as [START_REF] Lindsay | Inference functions and quadratic score tests[END_REF] refer the distance function in (2.6) as quadratic inference functions.

See also the comments on the GMM approach by Lindsay and Qu (2003, Section 2.3), for example.

Note that minimizing Q(β) in (2.6) with respect to β is equivalent to solve the GMM estimating equation

∂ψ 1 ∂β C 1 ψ 1 = 0, (2.8) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT for β, where ψ 1 = K -1 K i=1 ψ 1i with ∂ψ 1 ∂β = K -1 K i=1 T t=1 µ it x it x it .
As far as the C 1 matrix in (2.8) is concerned, by (2.7), this may be computed by

C -1 1 = K -2 K i=1 T u=1 T t=1 σ iut (β, σ 2 γ , ρ)x iu x it , (2.9)
where for all u, t = 1, . . . , T, σ iut (β, σ 2 γ , ρ) is given by (2.2) and (2.4).

GQL estimation for β

Denote the count response vector for the ith (i = 1, . . . , K) individual by y i = (y i1 , . . . , y it , . . . , y iT ) .

Let µ i (β, σ 2 γ ) = [µ i1 (β, σ 2 γ ), . . . , µ it (β, σ 2 γ ), . . . , µ iT (β, σ 2 γ )
] be the expectation of the response vector y i , where by (2.1)

µ it (β, σ 2 γ ) = exp(x it β + σ 2 γ /2). Next denote by Σ i (β, σ 2 γ , ρ) the covariance matrix of y i . Thus, Σ i (β, σ 2 γ , ρ) = (σ iut ), (2.10)
where the formulas for σ itt ≡ σ itt (β, σ 2 γ ) and σ iut ≡ σ iut (β, σ 2 γ , ρ) for u = t, are given by (2.2) and (2.4), respectively. Now, following [START_REF] Sutradhar | On Exact Quasilikelihood Inference in Generalized Linear Mixed Models[END_REF], for example, one may solve the generalized quasi-likelihood (GQL) estimating equation given by (2.11) [see [START_REF] Wedderburn | Quasilikelihood Functions, Generalized Linear Models, and the Gauss-Newton Method[END_REF] and [START_REF] Mccullagh | Quasilikelihood Functions[END_REF] for the independence case] to obtain the GQL estimate of β. For given σ 2 γ and ρ, the GQL estimate obtained from (2.11) is consistent for β. This is because, as E(Y i ) = µ i (β, σ 2 γ ), the estimating equation (2.11) is unbiased.

K i=1 ∂µ i (β, σ 2 γ ) ∂β Σ -1 i (β, σ 2 γ , ρ)[y i -µ i (β, σ 2 γ )] = 0,
In the next section, we compare the relative efficiency of the GQL estimate of β obtained from (2.11) to the GMM estimate obtained from (2.8).

Estimation of σ 2 γ

Notice in (2.2) that σ 2 γ plays a role for an over-dispersion index. Since in the mixed model set up, that is when ρ = 0 under the present set up, conditional on γ i , one can show that the count responses along with their squared and pairwise products are sufficient to estimate β and σ 2 γ (Jiang (1998), Sutradhar ( 2004)), we now exploit such information to estimate the σ 2 γ parameter. Note however that using the second order information in either GMM or GQL approach, will require the fourth order moments for the responses. Under the assumption that ρ = 0, these fourth order moments may be calculated directly. For example,

E(Y 2 iu Y iv Y it |ρ = 0) = E γ i [E(Y 2 iu |γ i )E(Y iv |γ i )E(Y it |γ i )] = E γ i [{µ * iu + µ * iu 2 }µ * iv µ * it ] = µ iu µ iv µ it exp(3σ 2 γ )[1 + µ iu exp(3σ 2 γ )].
(2.12)

GMM Estimation for σ 2 γ

In the fashion similar to that of the GMM estimation of β by solving (2.8), we obtain the GMM estimate for σ 2 γ by solving the GMM estimating equation

∂ψ 2 ∂σ 2 γ C 2 ψ 2 = 0, (2.13)
where

ψ 2 = K -1 K i=1 ψ 2i with ψ 2i = T u≤t [y iu y it -λ iut ], λ iut being the expectation of Y iu Y it for u ≤ t.
Note that by (2.2) and (2.4), we obtain

λ iut = E[Y iu Y it ] = µ it + [exp(σ 2 γ )]µ 2 it , for u = t ρ t-u µ iu + [exp(σ 2 γ )]µ iu µ it for u < t.
(2.14)

Further note that C 2 in (2.13) is not easy to compute for general ρ. We rather compute this matrix by using a 'working' longitudinal independence assumption, that is ρ = 0. The formula for C 2 is then given by

C -1 2 = var K -1 K i=1 ψ 2i = K -2 K i=1 T u≤ T m≤t φ iu mt (β, σ 2 γ , ρ = 0) , (2.15)
where for u ≤ and m ≤ t, the formulas for

φ iu mt (β, σ 2 γ , ρ = 0) = cov[Y iu Y it , Y im Y it |ρ = 0] = E[Y iu Y i Y im Y it |ρ = 0]-λ iu | ρ=0 λ imt | ρ=0
require the formulas for the fourth order moments that may be computed in the manner similar to that of (2.12).

GQL Estimation for σ 2 γ

Similar to the GMM estimation for σ 2 γ , in the GQL approach, we exploit the same basic data, namely the squared and pair-wise product of the observations, but in

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

a different way. Let the squared and the pair-wise product of the observations be denoted by u i = (y 2 i1 , . . . , y 2 iT , y i1 y i2 , . . . , y it y i,t+1 , . . . , y i,T -1 y iT ) .

By (2.14), the expectation of u i is given by

λ i (β, σ 2 γ , ρ) = E[U i ] = (λ i11 , . . . , λ itt , . . . , λ iT T , λ i12 , . . . , λ iut , . . . , λ i,T -1,T ) , (2.16)
where

λ itt ≡ λ itt (β, σ 2 γ ) = µ it + µ 2 it exp(σ 2 γ ), for t = 1, . . . , T, λ iut ≡ λ iut (β, σ 2 γ , ρ) = ρ t-u µ iu + µ iu µ it exp(σ 2 γ ), for u < t.
Now in the fashion similar to that of (2.11) for the GQL estimation of β, one may now write the GQL estimating equation for σ 2 γ as

K i=1 ∂λ i (β, σ 2 γ , ρ) ∂σ 2 γ Ω -1 i (β, σ 2 γ , ρ)[u i -λ i (β, σ 2 γ , ρ)] = 0,
(2.17) [START_REF] Sutradhar | On Exact Quasilikelihood Inference in Generalized Linear Mixed Models[END_REF]] where Ω i is the covariance matrix of u i . Note that it is however extremely cumbersome to compute the fourth order moments matrix, Ω i , under the auto-regression model (1.4). As a remedy, we use a 'working' covariance matrix of u i such as Ω iw (β, σ 2 γ , ρ = 0) for the general Ω iw (β, σ 2 γ , ρ), which would be computed by pretending that the repeated responses, conditional on γ i , are independent, even though in reality they are not. Thus, we propose to solve the 'working' GQL estimating equation

K i=1 ∂λ i (β, σ 2 γ , ρ) ∂σ 2 γ Ω -1 iw (β, σ 2 γ , ρ = 0)[u i -λ i (β, σ 2 , ρ)] = 0, (2.18)
for the estimation of σ 2 γ .

Estimation of the Longitudinal Correlation Parameter ρ

Note that the ρ parameter in the dynamic model (1.4) is understood in a much better way through the auto-covariances given in (2.4). It is clear from (2.4) that it would be sufficient to use the lag 1 pairwise product of the responses to estimate ρ. We exploit such information both in GMM and GQL approaches but in a different way.

The formulas for the estimate of ρ under these approaches are developed as follows.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

GMM Estimation for ρ

In the fashion similar to that of the GMM estimation of σ 2 γ by solving (2.13), we obtain the GMM estimate for ρ by solving the estimating equation

∂ψ 3 ∂ρ C 3 ψ 3 = 0, (2.19)
where

ψ 3 = K -1 K i=1 ψ 3i with ψ 3i = T t=2 [y it y i,t-1 -λ it,t-1 ],
where by (2.4)

λ it,t-1 = E[Y it Y i,t-1 ] = ρµ i,t-1 + [exp(σ 2 γ )]µ i,t-1 µ it .
(2.20)

Similar to the computation for C 2 in (2.15), we compute C 3 in (2.19) under the 'working' longitudinal independence assumption (ρ = 0). To be specific,

C -1 3 = var K -1 K i=1 ψ 3i = K -2 K i=1 T u=2 T t=2 φ iu,u-1,t,t-1 (β, σ 2 γ , ρ = 0) - K i=1 T t=2 λ it,t-1 (β, σ 2 γ , ρ = 0) 2   , (2.21)
where λ it,t-1 is given by (2.20), and the formula for φ iu,u-1,t,t-1 (β, σ 2 γ , ρ = 0) can be computed by using the formulas for the fourth order moments, similarly to (2.12).

ρ Estimation Under the GQL Approach

For given values of β and σ 2 γ , the correlation or probability parameter (ρ) may be consistently estimated by solving a suitable moment estimating equation that may be developed by equating the population covariance of the data given in (2.4), with its sample counterpart. Note that as ρ is a correlation parameter under the autoregressive order 1 set up, similar to the Gaussian set up, it would be sufficient to exploit the lag 1 auto-covariance only to estimate this parameter. More specifically, as by ( 2

.2), E(Y it -µ it ) 2 = σ itt = µ it +[exp(σ 2 γ )-1]µ 2 it , and by (2.4), E(Y it -µ it )(Y i,t+1 - µ i,t+1 ) = ρµ it + {exp(σ 2 γ ) -1}µ it µ i,t+1 , ρ may be estimated consistently by ρ = [a 1 -b 1 ]/g 1 , (2.22)
where a 1 is the observed lag 1 correlation defined as

a 1 = [ K i=1 T -1 t=1 y * it y * i(t+1) /K(T - 1)]/[ K i=1 T t=1 y * 2 it /KT ], with y * it = (y it -µ it )/(σ itt ) 1/2 , where σ itt = µ it + (exp(σ 2 γ ) - 1)µ 2 it . In (2.22), g 1 = K i=1 T -1 t=1 µ it (σ itt σ i,t+1,t+1 ) -1 2 /K(T -1), and b 1 = (exp(σ 2 γ ) - 1) K i=1 T -1 t=1 m it m i,t+1 /K(T -1), with m it = µ it /(σ itt ) 1 2 .
Since the means and the variances of the count responses are functions of both β and σ 2 γ , we are mainly interested to examine the relative efficiency of the GMM and GQL approaches in estimating these two parameters, by comparing their corresponding asymptotic variances.

Asymptotic variances of the GMM estimators

Since the K individuals provide repeated count responses independently, as K → ∞, it follows from (2.8) and (2.13), by using the central limit theorem that the GMM based estimators of β and σ 2 γ have the asymptotic covariances given by cov

( βGMM ) = lim K→∞ ∂ψ 1 ∂β C 1 ∂ψ 1 ∂β -1 (3.1) cov(σ 2 γ,GM M ) = lim K→∞ C -1 2 ∂ψ 2 ∂σ 2 γ -2 . (3.2)

Asymptotic variances of the GQL estimators

By similar calculations as in the GMM case, it follows from (2.11) that the asymptotic covariance matrix of the GQL estimator of β is given by cov

( βGQL ) = lim K→∞ K i=1 ∂µ i (β, σ 2 γ ) ∂β Σ -1 i (β, σ 2 γ , ρ) ∂µ i (β, σ 2 γ ) ∂β -1 . (3.3)
Next, it follows from (2.18) that the asymptotic variance of the GQL estimator of σ 2 γ is given by cov

(σ 2 γ,GQL ) = lim K→∞ K i=1 ∂λ i (β, σ 2 γ , ρ) ∂σ 2 γ Ω -1 iw (β, σ 2 γ , ρ = 0) ∂λ i (β, σ 2 γ , ρ) ∂σ 2 γ -1 . (3.4)

Asymptotic Efficiency Computation

We now illustrate the relative efficiency of the GQL and GMM estimators through a numerical example. For the purpose, we consider K = 500, p = 2, T = 4, and use a

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

covariate matrix with the first covariate as

x it1 =         
0 for i = 1, . . . , K/2; t = 1, 2 1 for i = 1, . . . , K/2; t = 3, 4 1 for i = K/2 + 1, . . . , K; t = 1, . . . , 4, . whereas the second covariate is chosen to be

x it2 =                1 for i = 1, . . . , K/2; t = 1, 2
1.5 for i = 1, . . . , K/2; t = 3, 4 0 for i = K/2 + 1, . . . , K; t = 1, 2 1 fori = K/2 + 1, . . . , K; t = 3, 4.

. Furthermore, for true parameter values, we consider β 1 = β 2 = 1.0; ρ = 0.3 and 0.8, and σ 2 γ = 0.5, 1.0, 1.5 and 2.0. The asymptotic variances computed from (3.1) through (3.4) are shown in Table 1 for the above selection of the parameter values.

[Insert Table 1 about here] It is clear from Table 1 that the variances of the estimators for all three main parameters β 1 , β 2 , and σ 2 γ , under the GQL approach are uniformly much smaller than the corresponding variances under the GMM approach, justifying that the GQL approach produces much more efficient estimates than the GMM approach for all main parameters of the model. For example, when ρ = 0.8 and σ 2 γ = 1.5, the GQL estimates of β 1 and β 2 are respectively 1.93×10.0 -2 6.62×10.0 -4 = 29.15 and 1.49×10.0 -2 4.86×10.0 -4 = 30.66 times more efficient than the corresponding GMM estimates. For the estimation of σ 2 γ , the GQL approach appears to perform extra-ordinarily better than the GMM approach.

For example, for the same set of parameters, i.e., when ρ = 0.8 and σ 2 γ = 1.5, the GQL estimate of σ 2 γ is 0.230 1.08×10.0 -5 = 21296.30 times more efficient. In summary, the GQL approach performs much better than the GMM approach in estimating all main parameters, its performance being extra-ordinarily better in estimating the variance component σ 2 γ .

ACKNOWLEDGMENT

This research was partially supported by a grant from the Natural Sciences and Engineering Research Council of Canada. The authors would like to thank a referee for comments and suggestions.

ACCEPTED MANUSCRIPT

Table 1. Comparison of asymptotic variances (Var) of the GQL and GMM estimators for the estimation of the regression parameters (β 1 and β 2 ), and the variance component (σ 2 γ ), of a longitudinal mixed model for panel count data, with T = 4 and K = 500.

Asymptotic variances ρ Method Quantity

σ 2 γ = 0.5 1.0 1.5 2.0 0.3 GQL Var( β1 ) 9.25 ×10.0 -4 7.88×10.0 -4 6.33×10.0 -4 4.99×10.0 -4 Var( β2 ) 7.11×10.0 -4 6.41×10.0 -4 5.27×10.0 -4 4.20 ×10.0 -4 Var(σ 2 γ ) 7.38 ×10.0 -5 3.03×10.0 -5 1.08×10.0 -5 1.04 ×10.0 -5 GMM Var( β1 ) 3.96 ×10.0 -3 9.63×10.0 -2 1.91×10.0 -2 3.47×10.0 -2 Var( β2 ) 2.99 ×10.0 -3 7.41 ×10.0 -3 1.47 ×10.0 -2 2.69×10.0 -2 Var(σ 2 γ ) 3.77 ×10.0 -3 3.08×10.0 -2 0.231 1.708 0.8 GQL Var( β1 ) 9.68 ×10.0 -4 8.21×10.0 -4 6.62×10.0 -4 5.24 ×10.0 -4 Var( β2 ) 6.74×10.0 -4 5.94×10.0 -4 4.86×10.0 -4 3.87×10.0 -4 Var(σ 2 γ ) 7.38×10.0 -5 3.03×10.0 -5 1.08×10.0 -5 1.04×10.0 -5 GMM Var( β1 ) 4.26×10.0 -3 9.86×10.0 -3 1.93×10.0 -2 3.48 ×10.0 -2 Var( β2 ) 3.16 ×10.0 -3 7.54×10.0 -3 1.49×10.0 -2 2.70×10.0 -2 Var(σ 2 γ ) 3.74×10.0 -3 3.07×10.0 -2 0.230 1.707


Derya Caliskan

Christophe Croux

Sarah Gelper

Efficient and robust scale estimation for trended time series

Keywords: Influence function, Robustness, Scale estimation, Time Series

This paper presents a new method for robust online variability extraction in time series. The proposed estimator is simultaneously highly robust and efficient. We derive its breakdown point, influence function, and asymptotic variance and study the finite sample properties in a simulation study.

Introduction

In the recent literature, new procedures have been proposed for robust scale estimation in a time series context, see for example [START_REF] Nunkesser | Online analysis of time series by the q n estimator[END_REF] and [START_REF] Gelper | Robust online scale estimation in time series: A regression-free approach[END_REF]. These procedures are designed to monitor the variability of noisy and trended time series. A drawback of the existing methods is that they lose efficiency with respect to the standard non-robust estimator. This paper proposes a new approach which combines the quantile based estimator of [START_REF] Gelper | Robust online scale estimation in time series: A regression-free approach[END_REF] and the efficient τ -scale estimator of

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT [START_REF] Yohai | High breakdown-point estimates of regression by means of the minimization of an efficient scale[END_REF]. As a result, we obtain a new estimator which is highly robust and at the same time attains high efficiency.

In several application fields data are automatically collected at a high frequency and need to be monitored instantaneously. Assume that we collect observations over time

y t = µ t + σ t t for t = 1, . . . , T (1)
where t is an i.i.d. sequence coming from a symmetric distribution F 0 with mean zero and variance 1. We estimate σ t using a moving window approach.

At each time t, a scale estimator σt is computed from a window of width n < T , containing the observations y t-n+1 , ..., y t . The variability is allowed to change slowly over time. More precisely, the variability is only assumed to be constant in a local window which can be fairly short, for example of only 20 observations.

For the ease of notation, we drop the time index t and the observations within a fixed window are denoted by y 1 , . . . , y n . We consider vertical heights of the triangles formed by triplets of successive data points. For three successive observations y i , y i+1 and y i+2 , the height of the triangle formed by these observations is given by [START_REF] Rousseeuw | Regression-free and robust estimation of scale for bivariate data[END_REF] propose robust scale estimators based on these heights in the context of nonparametric regression, and [START_REF] Gelper | Robust online scale estimation in time series: A regression-free approach[END_REF] adapted them to the time series context. These estimators are invariant if a trend line a + bt is added to the data. Moreover, they do not require to model and estimate the signal µ t , and were shown to be applicable for signals containing jumps, trends, trend changes and non-linearities.

h i = |y i+1 - y i + y i+2 2 | for i = 1, ..., n -2.
The estimator [START_REF] Gelper | Robust online scale estimation in time series: A regression-free approach[END_REF] advocate is the α-quantile of the heights of adjacent triangles

Q α adj (y 1 , ..., y n ) = c q h ( α(n-2)) , (2)
which is the α(n-2) -th value in the sequence of ordered heights, and where c q is a consistency factor. For α = 0.25 the highest value for the breakdown point is attained, and we denote Q adj ≡ Q 0.25 adj . While this estimator has good properties, it has a Gaussian efficiency of only 25%. [START_REF] Gelper | Robust online scale estimation in time series: A regression-free approach[END_REF] proposed to increase the value of α to get a better efficiency, at the price of a lower breakdown point. The procedure proposed in this paper, the τadjacent estimator, maintains the high breakdown point of the Q adj , but can attain an arbitrarily high efficiency.

This robust scale estimator τ adj is defined as

τ adj (y 1 , . . . , y n ) = c τ S 2 0 1 n -2 n-2 i=1 ρ( h i S 0 ) 1/2
(3) with initial scale estimator S 0 = Q adj (y 1 , . . . , y n ). The loss function ρ should be bounded, symmetric, ρ(0) = 0, and non-decreasing on the positive numbers. Its derivative ψ = ρ should exist and be positive at zero. In this paper we take Tukey's bisquare ρ function defined as

ρ(x) =    x 2 2 (1 -x 2 k 2 + x 4 3k 4 ), if |x| ≤ k k 2 6 if |x| > k.
The value of the tuning constant k depends on the desired level of efficiency.

As will be seen in Section 2, k = 5.48 yields a 95% asymptotic relative efficiency at the Gaussian model.

The performance of the adjacent τ -estimator is compared with the square root of the mean of squared adjacent heights,

M S adj (y 1 , . . . , y n ) = c s 1 n -2 n-2 i=1 h 2 i (4)
with c s again a consistency factor. The M S adj is not robust, but it is a standard proposal for scale estimation in nonparametric regression (see [START_REF] Gasser | Residual variance and residual pattern in nonlinear regression[END_REF]). It is not difficult to check that M S adj equals the τ adj for k → ∞.

Section 2 discusses statistical properties of the τ adj estimator. We show that the estimator has a breakdown point of 25%, a bounded influence function, and we compute its asymptotic variance. The good behavior of the estimator is confirmed in Section 3 by a simulation study. Section 4 summarizes the results.

Statistical Properties

In the section, we first derive the expression for the consistency factor c τ for the newly proposed estimator. Expressions for the constants c q in (2) and c s in (4) are given in [START_REF] Gelper | Robust online scale estimation in time series: A regression-free approach[END_REF]. We then derive the breakdown point and the influence function of the τ adj estimator. The breakdown point measures the robustness under larger amounts of outliers, while the influence function measures the sensitivity of the estimator with respect to small amounts of contamination. The last subsection computes the asymptotic variance of the estimator.

Fisher Consistency

For all theoretical derivations, we assume local linearity and a constant scale within the considered time window. Hence

y i = a + bi + σ i (5)
for i = 1, . . . , n where i ∼ F 0 . For an appropriately chosen window width

Let F be the distribution of the data and denote H F the distribution of the corresponding triangle heights. The functional form of the τ adj -scale

estimator (3) corresponds to τ adj (F ) = c τ S 2 0 (F )E H F ρ( h S 0 (F ) ) 1/2 (6)
where

S 0 (F ) ≡ Q adj (F ) = c q H -1 F (0.25).
The proposed scale estimator is location invariant, invariant when a trend is added to the data, and scale equivariant. Assume that the model ( 5) holds, and take a = b = 0 without loss of generality. Since the constant c q is such that Q adj is Fisher consistent, we have S 0 (F ) = σ. In order to achieve Fisher consistency, that is τ adj (F ) = σ, we need to take

c τ = 1 E H F 0 ρ( h S 0 (F 0 ) ) = 1 E H F 0 ρ(h) . ( 7
)
For F 0 a standard normal N (0, 1) distribution, it is not difficult to verify that E H F 0 ρ(h) = E F 0 ρ(Z 3/2), allowing for immediate calculation of (7).

As such, we obtain for k = 5.48 that c τ = 1.24. To make the estimator practically unbiased at finite samples, we propose to replace c τ by a finite sample version c n τ . By Monte Carlo simulation, following the approach outlined in [START_REF] Gelper | Robust online scale estimation in time series: A regression-free approach[END_REF], and for k = 5.48, we obtain

c n τ ≈ c τ n n -1.34
.

Breakdown point

The breakdown point of a scale estimator is the minimal amount of contamination such that the estimated scale becomes either infinite (explosion)

or zero (implosion). Let y n = {y 1 , . . . , y n } be a sample of size n. Let S be a scale estimator. Further, let y m n be a sample obtained from y n but with a proportion of m/n observations altered to arbitrary values (m ∈ {1, . . . , n}).

The finite sample breakdown point of S at the sample y n is defined as

ε * (S, y n ) = min 1 n m ∈ {1, 2, . . . , n} : sup y m n | log( S(y m n ) S(y n ) )| = ∞ ,
Suppose that y n is in general position, meaning that no three observations (i, y i ) lie on the same line for 1 ≤ i ≤ n. [START_REF] Gelper | Robust online scale estimation in time series: A regression-free approach[END_REF] have shown that, for the initial estimator

S 0 = Q adj , ε * (S 0 , y n ) = 1 n min n -1 -α(n -2) 3 , α(n -2) . (8)
The highest possible value for the breakdown point is attained for

α = n + 1 4(n -2) ≈ 0.25. ( 9
)
The finite sample breakdown point tends to the asymptotic breakdown point of 25%. In the following proposition we state that the breakdown point of the τ adj estimator is the same as that of the initial estimator S 0 .

Proposition 1. Let y n be a sample in general position. For the τ adj estimator with a bounded loss function ρ we have that

ε * (τ adj , y n ) = ε * (S 0 , y n ).

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Proof: From definition (3) of the τ adj estimator it is readily seen that the estimator tends to infinity if and only if the initial scale estimator S 0 does, since ρ is bounded. Implosion of the τ adj estimator occurs if either (i) S 0 implodes to zero (ii) all heights h i are equal to zero (iii) S 0 is so large that ρ(h i /S 0 ) is arbitrarily small for all 1 ≤ i ≤ n -2. Since we assumed that the sample is in general position, (ii) cannot occur. Furthermore, by definition of S 0 , about 75% of the heights is larger than S 0 , and ρ(1) > 0, hence also (iii) is excluded. We conclude that the τ adj implodes if and only of the initial scale estimator implodes.

Influence function

The influence function of the functional S at the distribution F measures the effect on S of adding a small mass at the point w, standardized by the mass of the contamination. If we denote the point mass at w by ∆ w and consider the contaminated distribution F ε = (1ε)F + ε∆ w then the influence function is given by

IF (w; S, F ) = lim ε→0 S(F ε ) -S(F ) ε . ( 10
)
In the Appendix we derive an explicit expression for the influence function of τ adj , assuming model ( 5) holds with F 0 = N (0, 1). Figure 1

Asymptotic variance

The estimators are based on heights of triangles. While the observations themselves are assumed to be independent, the heights will be autocorrelated up to order two. As in [START_REF] Genton | Asymptotic variance of M-estimators for dependent Gaussian random variables[END_REF], the asymptotic variance of an estimator S based on the heights h i is given by

ASV (S, F ) = E(IF 2 (h i ; S, H F )) + 2E(IF (h i ; S, H F )IF (h i+1 ; S, H F )) + 2E(IF (h i ; S, H F )IF (h i+2 ; S, H F )).
where IF (h; S, H F ) is the influence function of the estimator S at the distribution H F . At model ( 5) with F 0 = N (0, 1), the influence function of the heights can be obtained by straightforward calculus. Without loss of generality we may further assume that σ = 1. Then, )], Z ∼ N (0, 1), ψ = ρ , and N is an index referring to the assumption of normality. The influence function of the initial scale estimator Q adj is given by

IF (h; S, H N ) = IF (h; S 0 , H N ) c τ √ d - c τ 2 √ d E[ψ(Z)Z] + σ 0 2 c τ √ d ρ( h σ 0 ) - √ dc τ where σ 0 = 3 2 , d = E[ρ(Z
IF (h; S 0 , H N ) = c q [ 0.25 -I(h < Q α N ) 2 2/3ϕ( 2/3Q 0.25 N ) ] where Q 0.25 N = H -1 N (0.25)
is the first quartile of the distribution of the heights at the standard normal distribution, and I is the indicator function, see [START_REF] Gelper | Robust online scale estimation in time series: A regression-free approach[END_REF].

The exact value of the ASV for the non-trimmed mean-squared-heights estimator M S adj is 35/36, at the normal model with σ = 1. For the τ -scale estimator, the ASV is obtained by numerical integration. The asymptotic relative efficiency of the τ adj estimator w.r.t. M S adj is then defined as Eff(S, F ) = ASV (M S adj , F ) ASV (S, F ) .

Figure 2 present the asymptotic efficiency of the τ adj as a function of k, at the normal model. An efficiency of 0.95 is attained at k = 5.48. We also simulated, over M = 10000 simulation runs, finite-sample efficiencies for a window width of n = 20. One sees that the asymptotic results provide a good approximation for the finite sample setting. Furthermore, for k converging to zero, the efficiency of the Q adj and τ adj estimators coincide.

Simulation

In this section, we simulate a time series generated from model ( 1 where a patch is a groups of 3 consecutive outliers having the same value, and (c) innovation outliers. For schemes (a) and (b) we consider independent error terms ε t ∼ N (0, 1), and induce outliers by replacing a proportion ε of the observations, for ε = 0, 0.01, 0.05, and 0.10, by values coming from a N (0, 5). In simulation scheme (c) we consider a first order autoregressive model ε t = θε t-1 + v t , with v t ∼ N (0, 1), and θ = 0.5. Outliers are then induced by replacing a proportion ε of the v t by values coming from a N (0, 5).

For every simulated series, we compute the root mean squared error (RMSE)

RM SE = 1 T -n + 1 T t=n (σ t -σ t ) 2 σ 2 t 1/2 .
Here n = 20 is the window width, T is the length of the time series, and σt is the estimated scale.

In Table 1 the average RMSE over the M = 10000 simulated time series is reported. In absence of outliers, the robust τ adj is almost as efficient as

= c τ ∂ ∂ε [S 2 0 (F ε )E H F ε ρ( h S 0 (F ε ) )] 1/2 )| ε=0 = c τ 2 E H F ρ(h) -1/2 ∂ ∂ε [S 2 0 (F ε )E H Fε ρ( h S 0 (F ε ) )]| ε=0 = c τ 2 E H F ρ(h) -1/2     2IF (w; S 0 , F )E H F ρ(h) + ∂ ∂ε E H F ε ρ( h S 0 (F ε ) )| ε=0 A    
where we used that S(F 0 ) = 1. In the above expression, with ψ = ρ , (11)

A = ∂ ∂ε ∞ 0 ρ( h S 0 (F ε ) )dH F (h)| ε=0 = -IF (w; S 0 , F ) ∞ 0 ψ(h)hdH F (h) + ∞ 0 ρ(h)d ∂ ∂ε H F (h) B | ε=0 .
The above integrals can be computed either analytically or numerically. An expression for IF (w; Q adj , F ) is given in [START_REF] Gelper | Robust online scale estimation in time series: A regression-free approach[END_REF].

Figure 1 :

 1 Figure 1: Influence function of the τ adj estimator at the normal model with σ = 1.

Figure 2 :

 2 Figure 2: Asymptotic efficiencies (left) and finite sample efficiencies for n = 20 (right) of the τ adj estimator (solid line) and the Q adj estimator (dashed line) as a function of the tuning constant k at the normal model.

  τ adj , F ) = ∂ ∂ε τ adj (F ε )| ε=0

  One has H F (u) = Φ( 2/3u) -Φ(-2/3u), giving us for B:∂H F (u) ∂ε | ε=0 = -3(2Φ( 2/3h) -1) + Φ( √ 2(h + w)) -Φ( √ 2(wh))+ 2(Φ 4/5(w/2 + h)) -Φ 4/5(w/2h) := G(h, w).So we can writeA A = -IF (w; S 0 , F ) dG(h, w).We conclude that, withc 2 τ = 1/E H F 0 ρ(h), IF(w; τ adj , F ) = IF (w; S 0 , F ) 1 -

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Table 1: Simulated RMSE for clean data, 1%, 5% and 10% outliers, for window width n = 20, T = 1000, and averaged over M = 10000 simulation runs. Three types of outlier configurations are considered.

Additive outliers

Patches of outliers Innovation outliers the M S adj estimator. This is as expected, since the tuning constant k was selected to achieve a 95% relative efficiency. This in contrast with the Q adj estimator, having a RMSE which is almost twice as large. In presence of even only 1% of contamination, the M S adj is no longer the most precise, and the τ -estimator behaves best. For larger amounts of contamination, i.e. 10%, the Q adj estimator is slightly better than the τ adj , at least for additive outliers.

For patches of outliers and innovation outliers, the τ adj remains to have the smallest RMSE. We conclude from the simulation study that in presence of amounts of contamination up to 10%, the τ adj estimator is to be preferred.

If we have a larger (but less than 25%) proportion of outliers, then the bias of the τ adj remains bounded, given its high breakdown point, but the Q adj estimator tends to perform better. Since we expect that most univariate time series contain few outliers and that it is in practice quite rare to have large amounts of outliers, we recommend the τ adj estimator.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Conclusion

Robust scale estimators based on the heights of triangles formed by triplets of successive observation are used for monitoring the scale of nonlinear noisy time series, as is documented in [START_REF] Gelper | Robust online scale estimation in time series: A regression-free approach[END_REF]. In this paper we propose to use τ -estimators. These estimators keep the high breakdown point of the initial estimator, while they may have an arbitrarily high efficiency. The efficiency of the τ adj estimator depends on a tuning constant k. We computed the value of k yielding a 95% relative efficiency with respect to the standard estimator. Monte Carlo simulations illustrate the good performance of the proposed procedure.

A major question we did not addressed is the choice of the window width n. It needs to be small enough for (5) to hold, but large enough to still provide accurate estimates. Another topic for future research is to investigate the properties of the τ adj estimator for dependent data. The scale estimator can then still be applied, and will maintain the high breakdown point. The asymptotic variance, however, will depend on the dependency structure in the data.

Appendix

Derivation of the influence function of the τ adj estimator at the normal model:

Assume that model (5) holds, with F 0 = N (0, 1). Without loss of generality, assume a = b = 0, and σ = 1, such that F = N (0, 1). From (6) it follows


Kjell Konis
email: kjell.konis@epfl.ch

Konstantinos Fokianos
email: fokianos@ucy.ac.cy

Safe density ratio modeling

Keywords: biased sampling, differential expression, empirical likelihood, linear programming, separation

An important problem in logistic regression modeling is the existence of the maximum likelihood estimators.

Especially when the sample size is small, the maximum likelihood estimator of the regression parameters does not exist if the data are completely, or quasi-completely separated. Recognizing that this phenomenon has a serious impact on the fitting of the density ratio model-which is a semiparametric model whose profile empirical log-likelihood has the logistic form because of the equivalence between prospective and retrospective samplingwe suggest a linear programming methodology for examining whether the maximum likelihood estimators of the finite dimensional parameter vector of the model exist. It is shown that the methodology can be effectively utilized in the analysis of case control gene expression data by identifying cases where the density ratio model cannot be applied. It is demonstrated that naive application of the density ratio model yields to erroneous conclusions.

Consider two independent samples of observations such that X 11 , . . . , X 1n1 is a random sample from g 1 (x) = exp α + β T h(x) g 2 (x), X 21 , . . . , X 2n2 is a random sample from g 2 (x).

(1)

In the above g i (x), i = 1, 2 are unknown probability density functions, α is an unknown scalar parameter, β is a d-dimensional vector of parameters and h(x) is a d-dimensional vector which consists of known functions of X.

Model [START_REF] Konis | lpSolveAPI: Interface to lp solve 5[END_REF] is called density ratio model because it specifies that the log-ratio of two unknown probability density functions is linear in some parameters, see [START_REF] Anderson | Separate sample logistic discrimination[END_REF], [START_REF] Qin | A goodness of fit test for the logistic regression model based on case-control data[END_REF]. The model is motivated by means of the standard logistic regression and the equivalence between prospective and retrospective sampling, [START_REF] Prentice | Logistic disease incidence models and case-control studies[END_REF]. It should be mentioned that the class of distributions for which [START_REF] Konis | lpSolveAPI: Interface to lp solve 5[END_REF] holds is rather general and it includes the exponential family of distributions. Hence a vast collection of data types can be modeled by means of model ( 1).

An important observation is that when model ( 1) holds, and if β = 0, then the two samples are identically distributed. We conclude that model ( 1) is useful to the semiparametric comparison of two samples in the sense that the densities g i (.), i = 1, 2 are left completely unspecified but the weight function exp α + β T h(x) depends on some finite dimensional parameter. Hence (1) provides a sound framework for addressing the problem of comparing two independent samples. In addition, it provides a compromise between the fully parametric and nonparametric approaches to the problem of testing equality of two distribution, see [START_REF] Qin | Quantifying the change of melanoma incidence by Breslow thickness[END_REF], [START_REF] Kedem | Statistical comparisons of algorithms[END_REF] and [START_REF] Fokianos | A two-sample model for the comparison of radiation doses[END_REF], among others, for applications of the density ratio model to real data problems.

Inference regarding both finite and infinite dimensional parameters of model ( 1) has been studied by various authors assuming that the sample sizes tend to infinity in a suitable way. Following the empirical likelihood methodology, as advocated by [START_REF] Owen | Empirical Likelihood[END_REF], a parametric likelihood function for the finite dimensional parameters is obtained after profiling out the infinite dimensional parameter of the model. However, there are applications where the sample sizes are small and therefore direct application of the aforementioned techniques might suffer from non-existence or even non-convergence problems. To make this point clear, we follow [START_REF] Qin | A goodness of fit test for the logistic regression model based on case-control data[END_REF] and [START_REF] Fokianos | A semiparametric approach to the one-way layout[END_REF] who show that the empirical log-likelihood is given by

l(θ) ≡ l(α, β) = - 2 i=1 ni j=1 log 1 + ρ 1 exp(α + β T h(x ij )) + n1 j=1 α + β T h(x 1j ) (2)
where ρ 1 = n 1 /n 2 . Furthermore, if θ = (α, β) ′ denotes the maximum likelihood estimator of θ, assuming that it exists, then it can be shown that

pij = 1 n 2 1 1 + ρ 1 exp(α + βT h(x ij )) . ( 3
)
where p ij = dG 2 (x ij ), the size of the jump of G 2 (•) at the observed datum X ij = x ij . A consistent estimator for both of G 1 (.) and G 2 (.) can be constructed provided that the total sample size n = n 1 + n 2 tends to infinity such that n 1 /n 2 → ρ 1 -see [START_REF] Qin | A goodness of fit test for the logistic regression model based on case-control data[END_REF] and [START_REF] Fokianos | A semiparametric approach to the one-way layout[END_REF] for more.

Note that (2), after reparametrization, is equivalent to the standard logistic regression likelihood-a direct consequence of the equivalence between retrospective and prospective sampling as it was mentioned before. Hence, the finite dimensional vector of parameters θ = (α, β T ) T can be estimated by any of the numerous statistical programs which include logistic regression modeling. A standard approach for computing the maximum likelihood estimate of θ is to use the Fisher scoring method which, under some regularity assumptions, yields a sequence of approximations that converge to the maximum likelihood estimators of α and β. Occasionally this sequence of approximations does not converge to a finite value. Therefore application of the density ratio model is questionable, see [START_REF] Fokianos | Comparing two samples by penalized logistic regression[END_REF] who proposes penalization for the resolution of this problem and [START_REF] Davidov | On the existence and uniqueness of the npmle in biased sampling models[END_REF] who give necessary and sufficient conditions for the existence of maximum likelihood estimators for both the finite and infinite dimensional parameter. The non existence issue of the maximum likelihood estimators for the logistic regression model, sometimes, referred to as monotone likelihood or infinite parameters, occurs when a condition among the sample points, known as separation, prevails.

More specifically, the concept of separation in the logistic regression context was introduced by Albert and [START_REF] Albert | On the existence of maximum likelihood estimates in logistic regression models[END_REF] who showed that the sample points can be classified into one of three mutually exclusive configurations: complete separation, quasicomplete separation or overlap. To rephrase the concept of separation in terms of the density ratio model, define the sample points by the (d + 1)-dimensional vector

u ij = (1, h T (x ij )) T for j = 1, 2, . . . , n i , i = 1, 2.
Then, there exists complete separation among the sample points if there is a nonzero vector γ such that γ T u ij < 0 when i = 1 and γ T u ij > 0 when i = 2.

In other words, there is complete separation when there exits a hyperplane H such that all the values of the first sample lie strictly on one side of H and all the values of the second sample lie strictly on the other side of H. If complete separation is not present among the sample points but there exists a nonzero vector γ satisfying

γ T u ij ≤ 0 when i = 1 and γ T u ij ≥ 0 when i = 2,
then there is quasicomplete separation among the sample points. Finally, if there is no nonzero vector γ satisfying the last set of equations, then the sample points overlap. [START_REF] Albert | On the existence of maximum likelihood estimates in logistic regression models[END_REF] and [START_REF] Santner | Anderson's conditions for the existence of maximum likelihood estimates in logistic regression models[END_REF] show that there is a finite and unique maximum likelihood estimate of the logistic regression model parameter θ provided that there exists overlap among the sample points. Therefore, the density ratio model ( 1) is applicable only when there is overlap among the sample points according to the above conventions. The aim of this contribution is to suggest a test for checking whether there exists separation among sample points. By providing such a test, we answer the question of the applicability of model ( 1).

Methodology

The previous discussion shows that before fitting the density ratio model, we need to test for separation among the sample points. We outline a method based on linear programming. In the case that separation is found among the sample points we report that the density ratio model is not applicable. Otherwise (i.e., in the case of an overlap configuration), we fit the underlying logistic regression and proceed with inference, as in [START_REF] Fokianos | A semiparametric approach to the one-way layout[END_REF], for instance.

The use of linear programming to check for separated configurations of the sample points was first suggested in [START_REF] Albert | On the existence of maximum likelihood estimates in logistic regression models[END_REF] where the authors specified the necessary constraints but not an objective function.

Furthermore, [START_REF] Santner | Anderson's conditions for the existence of maximum likelihood estimates in logistic regression models[END_REF] described a mixed integer linear program capable of distinguishing between the three configurations while [START_REF] Silvapulle | Existence of maximum likelihood estimates in regression models for grouped and ungrouped data[END_REF] and [START_REF] Clarkson | Computing extended maximum likelihood estimates for linear parameter models[END_REF] used linear programming to check for the existence of a finite maximum likelihood estimate for the logistic regression model.

We test for separation based on recent work by [START_REF] Konis | Linear Programming Algorithms for Detecting Separated Data in Binary Logistic Regression Models[END_REF] because it is relatively simpler to apply. In addition, the test is implemented by means of the safeBinaryRegression R package and therefore it is easily accessible to data analysts, Konis (2009). In what follows, we assume that the n × (d + 1) design matrix U with rows given by u T ij is of full rank so that the associated logistic regression model is well-defined. Recalling the notation from the previous section, let s ij be the inner product of u ij and a vector γ = (γ 1 , . . . , γ d+1 ) T perpendicular to H, as shown in Figure 1. The linear program proposed in [START_REF] Konis | Linear Programming Algorithms for Detecting Separated Data in Binary Logistic Regression Models[END_REF] seeks to maximize the sum of the absolute values of the s ij subject to the constraint that s 1j ≤ 0 and s 2j ≥ 0. In the notation of linear (2001). These data were produced using cDNA technology and they describe intensity ratios of 3226 genes across 22 breast cancer samples. One of the aims of the experiment conducted by [START_REF] Hedenfalk | Gene-expression profiles in hereditary breast cancer[END_REF] was to detect differences in the genetic profiles between the BRCA1 and BRCA2 mutations. For these two conditions there are 7 and 8 arrays available, respectively. Notice that these data were log 10 transformed and therefore they assume both positive and negative values. We apply the density ratio model with h(x) = x (respectively, h(x) = x 2 ) and we obtain 19 (10, respectively) separated cases. In addition, the choice h(x) = (x, x 2 ) T yields 28 separated cases.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT h(x) x log x √ x (x, log x) T (x,
Some interesting observations follow next. Figure 3 shows data from a specific gene from the breast cancer study. For this particular case, it is worthwhile to notice that the density ratio model ( 1) is applicable when h(x) = x or h(x) = x 2 -a projection of the data onto either the x or the y axis produces an overlapped configuration among the sample points. However, the plot shows that when h(x) = (x, x 2 ) T then there is complete separation among the sample points and therefore model ( 1) is not applicable with this choice of h(•).

Figure 1 :

 1 Figure 1: Six sample points completely separated by the hyperplane H. The inner products s ij = u ij •γ are strictly positive when i = 2 and strictly negative when i = 1.

Figure 2 :

 2 Figure 2: Illustrating separation for seven genes in the Affymetrix Spike-In data. The filled triangles represent one sample and the up-side-down open triangles represent the other.

Furthermore, Figure 4 Figure 3 :

 43 Figure 3: The open circles correspond to gene 3094 of condition BRCA1 and the filled circles correspond to condition BRCA2 for the same gene. H is the separating hyperplane when h(x) = (x, x 2 ) T in (1).

Table 1 :

 1 x 2 ) T The number of separated cases for Affymetrix Spike-In study for different choices of h(•).

			# 15	15	15	556	1210
		36311_at
		36202_at
		36085_at
	Gene	33818_at
		32660_at
		1091_at
		1024_at
		1.9	2.0	2.1	2.2	2.3	2.4	2.5
					h(x) = log x

Acknowledgements

The authors thank A. Davison for constructive comments. This work was carried out while K. Fokianos was visiting the Institute of Mathematics, EPFL. The hospitality of all faculty members is greatly acknowledged.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

programming (see [START_REF] Konis | Linear Programming Algorithms for Detecting Separated Data in Binary Logistic Regression Models[END_REF]) this problem is stated by seeking the vector γ such that maximize: n2 j=1 u 2j • γ -n1 j=1 u 1j • γ, subject to: u 2j • γ ≥ 0, j = 1, 2, . . . , n 2 , u 1j • γ ≤ 0, j = 1, 2, . . . , n 1 , γ l ∈ R, l = 1, . . . , d + 1.

(4) If there is overlap among the sample points then there is no nonzero vector γ that satisfies the constraints imposed by (4). In this case γ ≡ 0 is the only feasible point and the optimal value of the objective function is zero. On the other hand, if there is either complete or quasicomplete separation among the sample points then there is γ such that at least one of the inner products is nonzero. However, if γ satisfies the constraints of (4) then so does kγ for real k > 1. Hence, in the separated case, the linear program is unbounded.

The safeBinaryRegression package uses lp solve (Berkelaar et al., 2009) via the lpSolveAPI (Konis and lp solve, 2009) R package to solve (4). If the linear program is bounded and the optimal value of the objective function is zero then there is overlap among the sample points and the density ratio model is appropriate. However, if the linear program is unbounded then the density ratio model cannot be applied to two sample problems.

Examples

Consider data from the Affymetrix Spike-In study [START_REF] Irizarry | Summaries of Affymetrix GeneChip probe level data[END_REF]. This experiment was used by Affymetrix to develop the MAS 5.0 preprocessing algorithms. The data set consists of measurements from 12626 human genes. We focus on the two array groups among the 14 array groups that contain 12 replicates each, leading to a case-control design with a total of 24 Affymetrix HGU95 biochips. For the data analysis, probe level summaries are computed using the RMA method [START_REF] Irizarry | Summaries of Affymetrix GeneChip probe level data[END_REF], yielding intensity levels on a log 2 scale. We test whether the density ratio model can be employed for inference for these data in the sense of comparing the distributions of cases and controls. Table 1 shows the number of separated cases for the whole data. It turns out that when h(•) is chosen to be univariate we obtain 15 separated samples and therefore the density ratio model is not applicable for these particular genes. The results are consistent for all choices of h(•) since both h(x) = log x and h(x) = √ x are monotone functions of x when x > 0. Note that the choices for h(•) are motivated by considering the log ratio of some standard distributions. For instance, choosing h(x) = x can be motivated by calculating the log ratio of two normal densities with identical variances. The non-applicability of the density ratio model turns out to be a more serious issue when h(•) is chosen to be a two dimensional function. In this case the number of separated samples increases considerably. Figure 2 demonstrates clearly the concept of separation for some selected genes from the Affymetrix Spike-In study.

As a second example we considered the data set described in the breast cancer study by Hedenfalk et al.

Conclusions

We have outlined a methodology for testing for separation among the sample points prior to fitting the density ratio model. It has been shown, using real data examples, that the proposed methodology effectively identifies cases where the model is not applicable. Some misleading results unfold if the technique is not applied properly.