A Chambaz

D Choudat

MAP5 C Huber

J-C Pairon

M J Van Der Laan

Threshold regression models adapted to case-control studies. Application to lung cancer induced by occupational exposure to asbestos

We rely on a recent French matched case-control study to investigate the effect of occupational exposure to asbestos on the occurrence of lung cancer. We build a large collection of threshold regression models, data-adaptively select a better model in it by multi-fold likelihood-based cross-validation, then fit the resulting better model by maximum likelihood. A necessary preliminary step to eliminate the bias due to the case-control sampling is made possible because the conditional distribution of being a case given the matching variable and the marginal distribution of the matching variable can be computed beforehand based on two studies independent from our dataset. The implications of the fitted model in terms of expected years of life free of lung cancer lost due to the occupational exposure to asbestos are discussed.

Introduction

Asbestos is a powerful carcinogen [IARC, 1977]. We rely on a recent French matched case-control study on lung cancer by [START_REF] Pairon | Occupational lung cancer: a multicentric case-control study in Paris area[END_REF] to investigate the effect of occupational exposure to asbestos on the occurrence of lung cancer. Following a case-control study design is convenient for a rare disease like lung cancer, with a known prevalence proportion approximately equal to five cases out of 10,000 persons [START_REF] Belot | Cancer incidence and mortality in France over the period 1980-2005[END_REF].

Logistic regression in parametric statistical models is the prevalent method of analysis in casecontrol study [Breslow, 1996, and references therein]. Sometimes however, as here due for instance to the nature of our exposure variable and because we do not aim for results in terms of risks, this approach would not be satisfactory. [START_REF] Gustavsson | Low-dose exposure to asbestos and lung cancer: Dose-response relations and interaction with smoking in a population-based case-referent study in stockholm, sweden[END_REF] recently undertook a casecontrol study on low-dose exposure to asbestos and lung cancer in Sweden with results expressed in terms of risks. By following van der Laan [2008], we show how to examine the connection between occupational exposure to asbestos and lung cancer, accounting notably for tobacco use, by performing estimation based on any parametric model for the conditional distribution of time to incident lung cancer given the remaining information. It is possible to take up that challenge because we can estimate the conditional distribution of being a case given the matching variable and the marginal distribution of the matching variable (based on the independent study undertaken by [START_REF] Belot | Cancer incidence and mortality in France over the period 1980-2005[END_REF] and on a dataset made publicly available by the French National Institute of Statistics and Economic studies, see http://www.insee.fr/en/).

We model the time to incidence of lung cancer as the time until an unobservable process crosses a threshold [START_REF] Lee | Threshold regression for survival analysis: modeling event times by a stochastic process reaching a boundary[END_REF]. The effect of occupational exposure to asbestos is included by accelerating the time index of the process. The effect of other covariates is not included as a time acceleration factor, but directly on the other process parameters. By separately modeling the effect of covariates and exposure in this way, we can get new insights on the effect of these covariates and also, more importantly, on the effect of exposure alone. We express the latter in terms of expected years of life lost as introduced by [START_REF] Robins | Estimability and estimation of expected years of life lost due to a hazardous exposure[END_REF]. Our model belongs to the family of threshold regression models, which has been playing an important role in survival analysis for some years [Lee andWhitmore, 2006, 2010, andreferences therein]. It is also an accelerated failure time model, with strong ties to the model developed by [START_REF] Oakes | Multiple time scales in survival analysis[END_REF] to study the effect of exposure to asbestos on the time until death from (not incidence of) lung cancer, as discussed later. We note that the study of a case-control dataset based on a threshold regression model is undertaken in [START_REF] Lee | A case-control study relating railroad worker mortality to diesel exhaust exposure using a threshold regression model[END_REF]], but we believe that the authors do not properly address the difficulties which stem from the type of sampling used to collect the data.

As our study involves an original qualitative description of the exposure to asbestos into 28 categories, we actually consider our model as a maximal model containing thousands of smaller models (obtained for instance by reducing the number of categories). We manage to data-adaptively select a better model in our large collection of threshold regression models by relying on multi-fold likelihood-based cross-validation [START_REF] Van Der Vaart | Oracle inequalities for multi-fold cross validation[END_REF]]. Then we fit the latter better model to the data by maximum likelihood, and draw our conclusions from its description.

The article is organized as follows. The dataset and the original qualitative description of the exposure to asbestos into 28 categories are described in Section 2. The case-control estimation problem is formalized in Section 3. We develop the threshold regression modeling in Section 4, showing how to derive the expected years of life free of lung cancer lost due to occupational exposure to asbestos for a case. We summarize in Section 5 the main results of the application to the dataset. A brief discussion is finally developed in Section 6. Some relevant material is gathered in the appendix. We present in Section A a few asymptotic results on the case-control weighted maximum likelihood estimator upon which our study relies. Elements of proofs are relegated to Section B. We provide additional details on the real data application in Section C. This includes the computation of the quantities required to eliminate the bias due to the casecontrol sampling in Section C.1, and brief presentations of (i) the principle of our multi-fold likelihood-based cross-validation model selection procedure in Section C.2, and (ii) its concrete implementation in Section C.3.

Dataset

In this section we describe the dataset of interest. The details of the case-control sampling scheme are given in Section 2.1. For each enrolled subject, non-professional information is collected, see Section 2.2, as well as an original history of occupational exposure to asbestos, see Section 2.3.

A matched case-control study

The matched case-control study took place between 1999 and 2002 in four Parisian hospitals. Case and control subjects were retrospectively recruited at the end of each year 1999 to 2002 among the patients of these hospitals who were free of lung cancer at the beginning of the corresponding year.

The case subjects were diagnosed with incident lung cancer during the period of the study. They were matched by control subjects on the basis of gender, age at end of calendar year (up to ±2.5 years), hospital, and race. Control subjects were recruited among patients of the departments of ophthalmology, general and orthopedic surgeries, and were by definition free of lung cancer at the time of their enrollment.

The one-to-one matching (i.e., the pattern of who is matched by whom) and race are not available. We propose an artificial valid matching pattern (based on gender, age and hospital) and make sure that our results do not depend on this particular choice. We exclude every subject with missing information. The resulting dataset counts n = 860 cases and 901 controls, hence a total of n + 901 = 1, 761 observations. We assume that the population sampled from during the study is stationary. Therefore, the observed data structures on experimental units made of pairs of case and matched control can be modeled as independent and identically distributed (iid) random variables. In Section 3, we invoke this fact to derive the likelihood function used in our study.

Non-professional information

Set a calendar time τ (expressed in years), and consider a subject sampled at time τ .

He/she is associated with his/her hospital of recruitment W 0 ∈ {1, 2, 3, 4}, gender W 1 (W 1 = 0 for men and W 1 = 1 for women), binary indicator W 2 of occurrence of lung cancer in close family (W 2 = 0 if no lung cancer occurred and W 2 = 1 otherwise), age at incident lung cancer diagnosis T (T = ∞ if no lung cancer occurred), and age at interview X = X(τ ).

It is well known that tobacco is a serious risk factor of lung cancer [START_REF] Biesalski | European consensus statement on lung cancer: risk factors and prevention. Lung cancer panel[END_REF]]. Thus information on tobacco consumption is collected during the interview. We summarize this information by considering a discretized version of the lifetime tobacco use. Hence, he/she is also associated with his/her tobacco use W 3 (W 3 = 0 for never-smoker, W 3 = 1 for lifetime tobacco use comprised between 1 and 25 pack years, W 3 = 2 for lifetime tobacco use comprised between 26 and 45 pack years, W 3 = 3 otherwise). The boundaries are chosen to yield strata of comparable sizes (371 subjects with W 3 = 0, and respectively 468, 469, 453 subjects with W 3 = 1, 2, 3). In particular, we overlook the duration of habit and years of abstinence, although they are known to play a role in the development of lung cancer [START_REF] Ruano-Ravina | Dose-response relationship between tobacco and lung cancer: new findings[END_REF][START_REF] Gustavsson | Low-dose exposure to asbestos and lung cancer: Dose-response relations and interaction with smoking in a population-based case-referent study in stockholm, sweden[END_REF].

We denote W = (W 1 , W 2 , W 3 ) ∈ W = {0, 1} 2 × {0, 1, 2, 3} the explanatory covariate, Z = min{T, X} and Y = 1{T ≤ X}. Note that Y = 1 if and only if (iff) T = Z ≤ X (the subject is then called a case) and Y = 0 iff T > Z = X (the subject is then called a control).

History of occupational exposure to asbestos

In addition, the occupational history up to age X, Ā(X), of the latter subject sampled at τ (we justify the notation below) is determined during the interview by experts of occupational exposure to asbestos. Every employment with duration at least 6 months is associated with its start and end dates, and with an original description of the exposure to asbestos. This description was characterized by the very experts who later conducted the interviews [START_REF] Pairon | Occupational lung cancer: a multicentric case-control study in Paris area[END_REF]. At the time of the characterization, the French threshold limit value for exposure to asbestos was set to 0.1f /mL (f /mL stands for "asbestos fibers per milliliter").

This description is a triplet referred to as "probability/frequency/intensity", each of them taking values in {1, 2, 3}: for the considered employment, the probability of exposure, its frequency and intensity are evaluated as either low or mild or high, respectively coded by 1, 2, 3. A probability index equal to 1, 2 or 3 corresponds to a passive exposure, a possible direct exposure or a very likely or certain direct exposure, respectively. A frequency index equal to 1, 2 or 3 corresponds to exposures occurring less than once a month, more than once a month and during less than half of the monthly working hours or during more than half of the monthly working hours, respectively. An intensity index equal to 1, 2 or 3 corresponds to a concentration of asbestos fibers less than 0.1 f/mL, between 0.1 and 1 f/mL and more than 1f/mL, respectively. Thus, the set E of categories of exposure has 27+1=28 elements (we add a category 0 for no exposure), each of them corresponding to a particular rate of exposure. Note that we will use from now on the notation ε

= ε 1 ε 2 ε 3 instead of ε = (ε 1 , ε 2 , ε 3 ).
The description is similar to that used by [START_REF] Gustavsson | Low-dose exposure to asbestos and lung cancer: Dose-response relations and interaction with smoking in a population-based case-referent study in stockholm, sweden[END_REF]. In the latter article, only the probability and intensity of exposure are considered. Four classes of probability of exposure correspond to an "estimated exposure prevalence" comprised between 0% and 20%, 20% and 50%, 50% and 85%, 85% and 100%. Four classes of intensity of exposure correspond to "time-period-specific annual arithmetic average level of exposure to asbestos" comprised between 0 and 0.03f/mL, 0.03f/mL and 0.1f/mL, 0.1f/mL and 0.3f/mL, 0.3f/mL and above. The Swedish threshold limit value for exposure to asbestos was set to 0.3f/mL in 1993.

We report in the top table of Table 1 the overall number of employments associated to each possible "probability/frequency/intensity" description. Although computed over a total of 8,432 employments, this table contains many zeros, showing that the latter description is over-parametrized. We also report in the bottom 1: Top: Overall number of employments associated to each possible "probability/frequency/intensity" description. The total number of employments is 8,432. Only 1,423 of them feature a description in E \ {0}. Bottom: Overall number of employments that feature a particular value of each coordinate of the "probability/frequency/intensity" description.

The generic longitudinal history of occupational exposure to asbestos is denoted by ā. It belongs to Ē, the set of functions from the nonnegative real line to E such that a(t) = 0 for t small or large enough (before the age at first employment or when no further information is available; this constraint is just a convenience, as we will make clear in Section 4). It is understood that the value of ā at t is denoted by a(t) while ā(t) stand for the restrictions of ā to [0, t]. Thus a(t) = a correspond to an occupational position held at age t and characterized by asbestos exposure a.

One of the central issues we deal with in this article is how to associate each description in E with a rate of exposure. We propose an original solution which heavily exploits the underlying multiplicative nature of the "probability/frequency/intensity" encoding. Indeed, it is the product of "probability", "frequency" and "intensity" which is relevant in terms of rate of exposure.

Formulation of the case-control estimation problem

We show in three steps how to examine the connection between occupational exposure to asbestos and lung cancer, accounting notably for tobacco use, by performing (appropriately weighted) maximum likelihood estimation based on any parametric model (i.e., in particular, not necessarily a logistic model) for the conditional distribution of time to incident lung cancer given the remaining information. First, we derive from the description of our case-control study the characterization of the representative sampling scheme one may have relied upon, had the probability of being an incident case of lung cancer not been so small. This mainly amounts to defining an observed data structure O ⋆ under representative sampling, whose distribution P ⋆ 0 presents features of interest, see Section 3.1. Second, we characterize the observed data structure O under matched case-control sampling and its distribution P 0 in terms of O ⋆ and P ⋆ 0 , see Section 3.2. Third, we show how to make inference on the features of P ⋆ 0 from data sampled under P 0 , see Section 3.3. This three-step procedure follows [van der Laan, 2008;[START_REF] Rose | Simple optimal weighting of cases and controls in case-control studies[END_REF].

Representative sampling scheme

We explained in Section 2.1 that sampling occurred at times τ 0 , τ 1 = τ 0 + 1, τ 2 = τ 0 + 2, τ 3 = τ 0 + 3 (where τ 0 stands for the initial sampling date, January 1st, 2000). Had the representative sampling been carried out, we would have observed n 0 (respectively n 1 , n 2 , n 3 ) independent observed data structures O ⋆ i sampled at time point τ 0 (respectively τ 1 , τ 2 , τ 3 ) under, say, P ⋆ (τ 0 ) (respectively P ⋆ (τ 1 ), P ⋆ (τ 2 ), P ⋆ (τ 3 )). We make the following stationarity assumption:

∀1 ≤ k ≤ 3, P ⋆ (τ k ) = P ⋆ (τ 0 ) ≡ P ⋆ 0 . ( 1
)
This assumption is justified by the influx and outflow experienced by the population of the Parisian region over the period of investigation. Thus, had the representative sampling been carried out, one would have finally collected

N = n 0 + n 1 + n 2 + n 3 independent copies (O ⋆ 1 , . . . , O ⋆ N ) of O ⋆ ∼ P ⋆ 0 , with O ⋆ = (W, X, Ā(X), Y, Z). (2)
Note that the likelihood of O ⋆ under P ⋆ 0 writes as

P ⋆ 0 (O ⋆ ) = P ⋆ 0 (W )P ⋆ 0 (X|W )P ⋆ 0 ( Ā(X)|X, W ) × dP ⋆ 0 (Z = T |T ≥ X -1, Ā(X), X, W ) Y × P ⋆ 0 (T > X|T ≥ X -1, Ā(X), X, W ) 1-Y , (3)
where dP ⋆ 0 (t|T ≥ X -1, Ā(X), X, W ) is the conditional density of T at time t given the event [T ≥ X -1, Ā(X), X, W ].

Matched case-control sampling

Such a representative sampling would have been impractical and ineffective because the probability P ⋆ 0 (Y = 1) of being an incident case of lung cancer is very small, see Section C.1. In order to recruit some cases in the sample, one would have to sample a huge number of observations. This is the main motivation for using a case-control sampling scheme.

We describe now what is our observed data structure in this framework. Introduce the categorical matching variable V ∈ V obtained by concatenating W 0 (subject's hospital when sampled), W 1 (subject's gender) and a discretized version of the age at sampling X over bins of length five years. In the sequel, we repeatedly use the convenient (though redundant) notation (V, W ).

The matched case-control sampling scheme can be described as follows:

• One first samples a case by sampling

(V 1 , O 1⋆ ) = (V 1 , W 1 , X 1 , Ā1 (X 1 ), Y 1 = 1, Z 1 )
from the conditional distribution of (V, O ⋆ ) given Y = 1 (the superscript "1" refers to the fact that Y 1 = 1).

• Subsequently, one samples J controls (V 0,j , O 0,j⋆ ) = (V 0,j , W 0,j , X 0,j , Ā0,j (X 0,j ), Y 0,j = 0, Z 0,j ) from the conditional distribution of (V, O ⋆ ) given Y = 0, V 0,j = V 1 for all j ≤ J (the superscript "0" refers to the fact that Y 0,j = 0 for all j ≤ J).

Conditional on V 1 = v ∈ V, the ratio of the number of controls for one case, J/1, is much smaller than the ratio

P ⋆ 0 (Y = 0|V = v)/P ⋆ 0 (Y = 1|V = v)
one would get in the population. This sampling scheme results in the observed data structure

O = ((V 1 , O 1⋆ ), (V 0,j , O 0,j⋆ ), j = 1, . . . , J) ∼ P 0
whose true distribution P 0 can be deduced from P ⋆ 0 and the two-step description above. The method naturally allows to consider the case that J is random and thus varies per experimental unit. This permits to exploit all our observations, even though we have less cases than controls. Note that each control is only taken into account once.

Case-control weighting of the log-likelihood loss function developed for representative sampling

It is remarkable that the log-likelihood loss function developed for the representative sampling scheme can be adapted to the case-control sampling scheme by appropriate weighting. This weighting relies on the prior knowledge of the joint distribution of (V, Y ), hence of the following probabilities: for each (y, v) ∈ {0, 1} × V,

q 0 = P ⋆ 0 (Y = 1), (4)
q 0 (y|v) = P ⋆ 0 (Y = y|V = v), (5) q 0 (v|y) = P ⋆ 0 (V = v|Y = y), (6)
or, namely, the marginal probability of being a case (4), the conditional probabilities of being a case or a control given matching variable at level v (5), and the conditional probabilities of observing level v for the matching variable given being a case or a control (6). We refer to Section C.1 for the computation of the latter key quantities. Consider a model P ⋆ for P ⋆ 0 . We assume without serious loss of generality that there exists a common dominating measure for all P ⋆ ∈ P ⋆ and P ⋆ 0 , so that we can refer to the densities p ⋆ and p ⋆ 0 of P ⋆ ∈ P ⋆ and P ⋆ 0 . With a slight abuse of notation, P ⋆ denotes both the model and the corresponding set of densities.

Define for all v ∈ V the quantities

q0 (v) = q 0 P ⋆ 0 (Y = 0|V = v) P ⋆ 0 (Y = 1|V = v) = q 0 q 0 (0|v) q 0 (1|v) , (7)
and introduce the following case-control weighted log-likelihood loss function for the density p ⋆ 0 of P ⋆ 0 under sampling of O ∼ P 0 : for all p ⋆ ∈ P ⋆ ,

ℓ(p ⋆ )(O) = q 0 log p ⋆ (V 1 , O 1⋆ ) + q0 (V 1 ) 1 J J j=1 log p ⋆ (V 1 , O 0,j⋆ ). ( 8
)
Even though q 0 appears in both terms in (8), we prefer to consider ℓ(p ⋆ )(O) as defined above rather than q -1 0 ℓ(p ⋆ )(O). This choice guarantees that the weighted log-likelihood ℓ(p ⋆ )(O) is on the same scale as the log-likelihood log P ⋆ 0 (O ⋆ ) under representative sampling. The weighted loglikelihood loss function is adapted to the case-control sampling scheme in the following sense (the proof is relegated to Section B):

Proposition 1. Let us assume that model P ⋆ for p ⋆ 0 is such that log p ⋆ (o ⋆ )dP ⋆ 0 (o ⋆ , Y = y)
are properly defined for all p ⋆ ∈ P ⋆ and y = 0, 1. If P ⋆ is well-specified (i.e., if p ⋆ 0 ∈ P ⋆ ), then the density that maximizes the expectation under P 0 of the weighted loss function (8) over P ⋆ , arg max

p ⋆ ∈P ⋆ E P0 ℓ(p ⋆ )(O),
is unique and coincides with p ⋆ 0 . In words, this tells us that it is possible to draw inference about P ⋆ 0 based on data sampled from P 0 by (appropriately weighted) maximum likelihood.

In the next section, we propose a parametric model for the conditional distribution of O ⋆ given (W, X, Ā, Y ), that is the conditional distribution of Z given (W, X, Ā, Y ). The parametric model is sound in the sense that the conditional distribution of Z given (W, X, Ā, Y ) only depends on Ω = (W, X, Ā(X), Y ). The latter parametric model is combined with a nonparametric model for the conditional distribution of (W, X, Ā) given Y , both yielding a semiparametric model P ⋆ = {P ⋆ θ : θ ∈ Θ} for P ⋆ 0 because we know beforehand q 0 , the true probability of being a case. Formally, the likelihood of O ⋆ = (Ω, Z) under θ ∈ Θ writes as

P ⋆ θ (O ⋆ ) = p ⋆ θ (Z|Ω)η(Ω),
η(Ω) being the likelihood of Ω, which we assume without serious loss of generality to be bounded away from 0. If we set

Ω 1 = (W 1 , X 1 , Ā1 (X 1 ), Y 1 = 1), Ω 0,j = (W 0,j , X 0,j , Ā0,j (X 0,j ), Y 0,j = 0), hence O 1⋆ = (Ω 1 , Z 1
) and O 0,j⋆ = (Ω 0,j , Z 0,j⋆ ), then the weighted log-likelihood loss function for p ⋆ 0 under sampling of O ∼ P 0 satisfies, for all θ ∈ Θ,

ℓ(p ⋆ θ )(O) = q 0 log p ⋆ θ (V 1 , O 1⋆ ) + q0 (V 1 ) 1 J J j=1 log p ⋆ θ (V 1 , O 0,j⋆ ) = q 0 log p ⋆ θ (Z 1 |Ω 1 , V 1 ) + q0 (V 1 ) 1 J J j=1 log p ⋆ θ (Z 0,j |Ω 0,j , V 1 ) +rem(O), (9)
where rem(O) is a random term independent of θ. We set l(θ

)(O) = ℓ(p ⋆ θ )(O) -rem(O) for all θ ∈ Θ.
Assuming that the true density p ⋆ 0 is "projected" (in terms of Kullback-Leibler divergence) onto p ⋆ θ0 for some θ 0 ∈ int(Θ), or in other terms that the mapping θ → KL(p ⋆ 0 , p ⋆ θ ) achieves a unique minimum at the unique θ 0 ∈ int(Θ), we focus hereafter on the maximum likelihood estimation of θ 0 . Following the lines of the proof of Proposition 1, the case-control weighted maximum likelihood estimator

θ n = arg max θ∈Θ n i=1 l(θ)(O i ) (10)
does estimate θ 0 . We briefly consider some asymptotic properties (including consistency and asymptotic normality) in Section A.

Threshold regression parametric modeling

We model the time to incidence of lung cancer as the time until an unobservable process crosses a threshold, the effect of asbestos exposure being included by accelerating the time index of the process, see Section 4.1. This is typical of threshold regression modeling [Lee andWhitmore, 2006, 2010, andreferences therein], with an original acceleration tailored to the description of occupational exposure to asbestos (as presented in Section 2.3) that we introduce in Section 4.2. We complete the description of our model in Section 4.3, where we also specify the form of the resulting case-control weighted log-likelihood loss function. Finally, we show in Section 4.4 how our model yields insights on the association between asbestos exposure alone and risk of developing lung cancer.

Time to incident lung cancer

Let B be a Brownian motion. For any real numbers h > 0 and µ ≤ 0, define

T {h, µ} = inf{t ≥ 0 : h + µt + B t ≤ 0}, ( 11
)
the first time the drifted Brownian motion (h + µt + B t , t ≥ 0) crosses the threshold 0. The distribution of T {h, µ} is known as the inverse Gaussian distribution with parameter (h, µ). It is characterized by its cumulative distribution function (cdf): for all t > 0,

F (h, µ)(t) = 1 + e -2hµ Φ (µt -h)t -1/2 -Φ (µt + h)t -1/2 ,
where Φ is the standard normal cdf. It is well known [START_REF] Chhikara | The inverse Gaussian distribution: theory, methods and applications[END_REF]] that T {h, µ} < ∞ almost surely eventually because µ ≤ 0. Therefore the distribution of T {h, µ} is also characterized by its density

f (h, µ)(t) = h (2πt 3 ) 1/2 exp - (h -|µ|t) 2 2t (all t > 0). Finally, T {h, µ} has mean h/|µ| whenever µ < 0.
Here, (h + µt + B t , t ≥ 0) models the amount of health relative to lung cancer in absence of occupational exposure to asbestos. Describing what happens in presence of exposure to asbestos involves the introduction of an acceleration function R which summarizes the effects of this exposure. An acceleration function R is a nondecreasing continuous function on the nonnegative real line such that R(t) ≥ t for all t ≥ 0. Given such a function R, we define

T {h, µ, R} = inf{t ≥ 0 : h + µR(t) + B R(t) ≤ 0}, ( 12
)
the first time the drifted Brownian motion (h + µt + B t , t ≥ 0) crosses the threshold 0 along the modified time scale derived from R. Obviously, T {h, µ, R} = T {h, µ} when R is the identity, but in general T {h, µ, R} ≤ T {h, µ}. Furthermore,

T {h, µ, R} ≥ t iff T {h, µ} ≥ R(t), (13)
so that the cdf of T {h, µ, R} at t > 0 is F (h, µ)(R(t)), and its density at

t > 0 is R ′ (t)f (h, µ)(R(t))
as soon as R is differentiable. Consequently, the conditional survival function and density of

T {h, µ, R} at t ≥ x -1 given [T {h, µ, R} ≥ x -1] are respectively G(h, µ, R)(t) = 1 -F (h, µ)(R(t)) 1 -F (h, µ)(R(x -1)) (14) and g(h, µ, R)(t) = R ′ (t)f (h, µ)(R(t)) 1 -F (h, µ)(R(x -1)) , (15)
two important equalities in view of the factorization of the likelihood shown in (3). Parameters h, µ and T {h, µ} are respectively interpreted as an initial amount of health relative to lung cancer, a rate of decay of this amount of health in absence of occupational exposure to asbestos, and the time to incident lung cancer in absence of occupational exposure to asbestos. As for T {h, µ, R}, it is interpreted as the time to incident lung cancer for a history of occupational exposure to asbestos summarized by acceleration function R. Admitting that the reference time scale (that is that of the Brownian motion B) corresponds to chronological/calendar time scale, the new time scale formed by the acceleration function R may be understood as a biological time scale. This interpretation acknowledges the fact that the ageing phenomenon related to lung cancer is stronger in presence of noxious occupational exposure to asbestos than in its absence.

Calendar versus biological ages: modeling the ageing acceleration due to occupational exposure to asbestos

We present now an original class of acceleration functions tailored to our particular description of occupational exposures to asbestos. Let us define

M = (M 0 , (M k,l ) k,l≤3 ) ∈ R + × M 3,3 (R + ) : 0 ≤ M k,1 ≤ M k,2 ≤ M k,3 = 1, k = 1, 2, 3 . (16)
Then the rate yielded by description

ε = ε 1 ε 2 ε 3 ∈ E \ {0} for M ∈ M writes as M (ε) = 1 + M 0 × M 1,ε1 × M 2,ε2 × M 3,ε3 , (17)
and that of ε = 0 is set to M (0) = 1. Notably, M 0 is the factor of acceleration of time for the higher exposure, which we recall is encoded by ε = 333. Rates M (ε) range from 1 to M (333) = 1 + M 0 and (with convention 0/0 = 1)

M (ε) -1 M 0 = M 1,ε1 × M 2,ε2 × M 3,ε3 :
an exposure characterized by "probability/frequency/intensity" description

ε = ε 1 ε 2 ε 3 achieves a fraction M 1,ε1 × M 2,ε2 × M 3,ε3
of the maximal acceleration. Note that we only need 7 parameters in order to fully describe the 28 possibly different rates of acceleration. Furthermore, it is easily seen that this parametrization is identifiable

: if M, M ′ ∈ M satisfy M (ε) = M ′ (ε) for all ε ∈ E then M = M ′ .
Consider M ∈ M and a generic longitudinal history ε as presented in Section 2. The mapping M (ε(•)) : t → M (ε(t)) is piecewise constant but it is more convenient to consider a continuous approximation to it, which we denote by r(M, ε). (Formally, one may rely on convolution to derive r(M, ε) from M (ε(•)).) Finally, every pair (M, ε) gives rise to the acceleration function R(M, ε) characterized by

R(M, ε)(t) = t 0 r(M, ε)(s)ds (18)
(all t ≥ 0). The quantity R(M, ε)(t) is a summary measure of the history of occupational exposure to asbestos, and can be interpreted as the cumulative occupational exposure to asbestos up to age t. In particular if ε(t) = 0 for all t ≥ 0 (i.e., in absence of occupational exposure to asbestos throughout lifetime) or if M 0 = 0 (i.e., assuming that exposure to asbestos has no effect on the risk of developing lung cancer), then R(M, ε)(t) = t for all t ≥ 0: in other words, the chronological and biological time scales coincide. Note that

R(M, ε) is differentiable because r(M, ε) is continuous. Note that R(M, ε)(t)
is very close to a linear combination of the times spent in each job category (if larger than 6 months) up to calendar time t, where the coefficients of the combination depend through M on the job categories. In this view, our acceleration function is classical [START_REF] Lee | Threshold regression for survival analysis: modeling event times by a stochastic process reaching a boundary[END_REF]. What makes it original though is how M maps each description ε ∈ E to its coefficient M (ε). We comment further on our acceleration function in Section 4.3. Now, given parameters h > 0, µ, M ∈ M and covariate ā, we obtain R a = R(M, ā), which yields in turn the time to incident lung cancer T {h, µ, R a } for the history of occupational exposure to asbestos summarized by R a .

The parametric model, and related case-control weighted log-likelihood loss function

For this purpose, we state that under θ = (α, β, M ) ∈ Θ = R 4 × R 16 × M, the conditional distribution of T (the possibly unobserved time to incident lung cancer of the subject associated with

O ⋆ ) given Ω is that of T {h, µ, R a } with log h = α index(W1,W2) (19)
(each level of (W 1 , W 2 ) is associated with a unique positive initial health h),

log(-µ) = β index(W1,W2,W3) (20)
(each level of (W 1 , W 2 , W 3 ) is associated with a unique negative drift µ), and

R a = R(M, Ā(X)).
Therefore, it holds that log 14) and ( 15) for the definitions of G and g). Finally, the relevant part of the resulting case-control weighted log-likelihood at θ ∈ Θ writes as

p ⋆ θ (Z|Ω) = Y log g(θ)(Z) + (1 -Y ) log G(θ)(Z), with convention G(θ)(Z) = G(h, µ, R a )(Z) and g(θ)(Z) = g(h, µ, R a )(Z) (see (
n i=1    q 0 log g(θ)(Z 1 i ) + q0 (V 1 i ) 1 J i Ji j=1 log G(θ)(Z 0,j i )    = P n l(θ), (21)
where l(θ 9)) and P n = n i=1 δ Oi denotes the empirical measure.

)(O) = ℓ(p ⋆ θ )(O) -rem(O) (see equation (
Comments. [START_REF] Oakes | Multiple time scales in survival analysis[END_REF] carried out the study of the effect of exposure to asbestos on time until death from (not incidence of) lung cancer. Although the path that leads us to our model is quite different from his (asbestos exposure are measured differently, more covariates are accounted for here than there), the conditional cdf at t ≥ 0 of time until the event of interest given the rest ends up in both cases being of the form F (R a (t)). Here F = F (h, µ) and R a = R(M, Ā(X)) whereas there F is a Weibull cdf and R a : t → t + ρcum(t), cum(t) representing the cumulative exposure to asbestos at age t in 100 million particles per cubic foot (measured in mppcf-years), and ρ an "equivalence parameter". Hence both model are accelerated failure time models of the simple collapsible form (since F (R a (t)) depends on the cumulated exposure at time t and not on the entire exposure history up to t; see [START_REF] Duchesne | On the collapsibility of lifetime regression models[END_REF] for theoretical conditions under which an accelerated failure time model turns out to be collapsible). Note that the expression "equivalence parameter" conveys the notion that a cumulative exposure of 1/ρ mppcf-years of asbestos dust has the same effect on the cumulative lung cancer risk as an additional year of life. In view of ( 17) and ( 18), the two acceleration functions are very similar, the component

M 0 of M = (M 0 , (M k,l ) k,l≤3
) playing the same role as ρ, and (M k,l ) k,l≤3 serving as a mean to associate a concentration of exposure to asbestos to every ε ∈ E.

Furthermore, we emphasize that our case-control weighted log-likelihood differs from the loglikelihood developed in [START_REF] Lee | A case-control study relating railroad worker mortality to diesel exhaust exposure using a threshold regression model[END_REF] even though the latter article also relies on a case-control study. [START_REF] Lee | A case-control study relating railroad worker mortality to diesel exhaust exposure using a threshold regression model[END_REF] justify the form of their log-likelihood and discuss on the limitations of the inferences they subsequently draw from it. However, we believe that the difficulties that stem from the fact that a case-control sampling is performed are not fully addressed by them.

Expected years of life free of lung cancer lost due to occupational exposure to asbestos

Equivalence (13) has an important consequence: given parameters h > 0, µ ≤ 0, M ∈ M and history of occupational exposure to asbestos Ā(X), T {h, µ} = R a (T {h, µ, R a }). In words, all things (gender, occurrence of lung cancer in close family, lifetime tobacco use) being equal, the time to incident lung cancer in the absence of occupational exposure to asbestos can be deduced deterministically from the (observed) age at incident lung cancer and history of occupational exposure to asbestos of a case. The nonnegative quantity

R a (T {h, µ, R a }) -T {h, µ, R a } (with convention R a (∞) = ∞ and ∞ -∞ = 0
) can be interpreted as a number of years of life free of lung cancer that the case could have enjoyed had he/she not been exposed to asbestos. Note that R a (T {h, µ, R a }) -T {h, µ, R a } is different from the remaining number of years of life free of lung cancer, as death may occur anytime after T {h, µ, R a } even in the absence of occupational exposure to asbestos.

There is a strong connection between R a (T {h, µ, R a }) -T {h, µ, R a } and expected years of life lost as introduced by [START_REF] Robins | Estimability and estimation of expected years of life lost due to a hazardous exposure[END_REF]. Following [START_REF] Robins | Estimability and estimation of expected years of life lost due to a hazardous exposure[END_REF], consider the expected years of life lost ∆(t) for a randomly sampled person among subjects (i) who share the same characteristics (gender, occurrence of lung cancer in close family, lifetime tobacco use, hence the common parameter (h, µ)), (ii) who develop an incident lung cancer at age t, and (iii) whose histories of occupational exposure to asbestos up to age t coincide (hence the same acceleration function R a , at least up to time t). Then under some assumptions on how the exposure affects health [Robins and Greenland, 1991, Theorem 3], it holds that ∆(t) = S -1 0 •S 1 (t)-t, where

S 0 = 1 -F (h, µ) and S 1 = (1 -F (h, µ)) • R a , so that ∆(t) = R a (t) -t.
For this reason, we decide to refer to R a (T {h, µ, R a }) -T {h, µ, R a } as expected years of life (free of lung cancer) lost (due to occupational exposure to asbestos).

Results

We present here the inferences we draw from our dataset, model, and resulting case-control weighted log-likelihood loss function as described in Sections 2 and 4.3. In Section 5.1 we comment on the fitted model, then we focus in Section 5.2 on its implications in terms of expected years of life lost.

Fitting the best model

It makes no doubt that the model we have built so far is over-dimensioned. The "probability/frequency/intensity" description with its 28 different levels is itself certainly too rich (see Table 1), or at least difficult to establish and prone to errors. We rather consider the model {P ⋆ θ : θ ∈ Θ} described so far as a "maximal" model giving rise to a large collection of sub-models {P ⋆ θ : θ ∈ Θ k } obtained by adding constraints on the "maximal" parameter θ = (α, β, M ) ∈ Θ. The number of such sub-models is large indeed: there are (1 + 7 3 ) = 344 sub-models defined by adding only constraints on M (of the type M 0 = 0, or M 0 > 0 and for any k

= 1, 2, 3, 0 = M k,1 or M k,1 = M k,2 or M k,2 = 1 or 0 = M k,1 = M k,2 or M k,1 = M k,2 = 1 or (0 = M k,1 , M k,2 = 1)
), hence the total number of sub-models equals 2 2 × 2 3 × 344 = 11, 008. It is out of question to explore the whole collection of sub-models. Instead, we propose to (i) define a large collection {Θ k : k ∈ K} of interest, (ii) let the data select a better Θ k in the latter collection based on a multi-fold likelihood-based cross-validation criterion.

van der Vaart et al. [2006] show that, under mild assumptions, the multi-fold likelihood-based cross-validation criterion will select a better model comparing favorably with the oracle model of the collection (whose definition involves the true distribution of the data). By this we mean that the likelihood risk of the better model will not be much bigger than that of the oracle model. Although we cannot invoke rigorously this remarkable property here, it motivates the procedure that we describe in Sections C.2 and C.3. The best model {P ⋆ θ : θ ∈ Θ k} is described in Section C.3. Its characterization teaches us that neither the initial health nor the drift parameter depend on the indicator of occurrence of lung cancer in close family. Moreover, (i) a passive exposure (probability index equal to 1) is the same as no exposure at all, and (ii) being exposed less than once a month (frequency index equal to 1) or between once a month and during less than half of the monthly working hours (frequency index equal to 2) have the same effect.

We first fit the best model in terms of maximum likelihood on the whole dataset. Regarding the derivation of confidence intervals, we decide to rely on the bootstrap instead of a central limit theorem (such as Proposition 3 in Section A). The particulars of the bootstrap procedure follow. We set α = 2.5%, B = 1, 000 and p = 5%, then for b ranging from 1 to B, we repeatedly resample without replacement n(1p) = 817 observed data structures, yielding the bootstrapped empirical measure P b n(1-p) , in order to compute and store the corresponding maximum likelihood estimate k are used as lower-and upper-bounds of confidence intervals, which simultaneously provide a (1 -2α) = 95%-coverage by the applied Bonferroni correction. Specifically, we obtain:

θ n(1-p), k(P b n(1-p) ) of θ ∈ Θ k. The
• initial health: 82 [23.42; 24.13] 1 25. 09 [24.86; 25.40] It is seen in particular that women are associated with a significantly larger initial health than men.

W 1 h 0 23.
• drift: 89 [13.25; 14.46] 10.55 [9.63; 11.80] 3 17.67 [17.11; 18.38] 14.79 [13.65; 17.77] Two main features arise:

-100µ W 3 W 1 = 0 W 1 = 1 0 0.
-For each level of lifetime tobacco use, the absolute value of the drift is significantly larger for men than for women (actually, the confidence intervals for W 3 = 3 slightly overlap). Combined with the already mentioned fact that women are associated with a larger initial health, this implies that for any given history of exposure to asbestos and for every level of lifetime tobacco use, the distribution of time to incident lung cancer in women is stochastically dominated by the distribution of time to incident lung cancer in men. In other words, given a man and a woman sharing the same history of exposure to asbestos and lifetime tobacco use, given an age t, the man is more likely to have developed an incident lung cancer at age t than the woman. Note that there is no clear consensus in the literature on whether there exist differences in lung cancer risk between men and women or not (for instance, [START_REF] Zang | Differences in lung cancer risk between men and women: examination of the evidence[END_REF] argue that women are more susceptible to tobacco carcinogens, but [START_REF] Haiman | Ethnic and racial differences in the smoking-related risk of lung cancer[END_REF] show that men or women are more susceptible to tobacco carcinogens, depending on ethnic and racial group).

-Both in men and women, the absolute value of the drift significantly increases with lifetime tobacco use. This implies that, both in men and women, for any given history of exposure to asbestos and for every 0 ≤ w < w ′ ≤ 3, the distribution of time to incident lung cancer for lifetime tobacco use equal to w is stochastically dominated by .192 [0.344; 1.998] Table 2: Estimated values (precision 10 -3 ) of the factor of acceleration of time (M (ε) -1) and related confidence intervals for each level of exposure ε ∈ E \ {0}. Recall that M (0) = 1.

ε M (ε) -1 ε M (ε) -1 ε M (ε)
the distribution of time to incident lung cancer for lifetime tobacco use equal to w ′ . In other words, given two persons sharing the same gender and history of exposure to asbestos, the person with the larger lifetime tobacco use is more likely to have developed an incident lung cancer at age t than the other. This is in agreement with the general scientific consensus [START_REF] Biesalski | European consensus statement on lung cancer: risk factors and prevention. Lung cancer panel[END_REF]].

• exposure to asbestos:

M 0 : 1.19 [0.34; 2.00] M 1,1 = 0 M 1,2 : 0.97 [0.96; 0.99] M 1,3 = 1 M 2,1 = M 2,2
M 2,2 : 0.93 [0.90; 0.98] M 2,3 = 1 M 3,1 : 0.02 [0.00; 0.09] M 3,2 : 0.09 [0.00; 0.27] M 3,3 = 1 We notably derive from the above table the values of (M (ε)-1) (which can be interpreted as a factor of acceleration of time due to an exposure of level ε, see ( 17)) and related confidence intervals for each level of exposure ε ∈ E \ {0}, see Table 2.

Application to the expected years of life free of lung cancer lost due to occupational exposure to asbestos

In view of Section 4.4, the results of the previous section provide us with a way of evaluating the expected years of life (free of lung cancer) lost (due to occupational exposure to asbestos) on a case by case basis. Say that we mostly care for a pointwise estimation of, and confidence lowerbound on, the expected years of life lost. A confidence upper-bound could be derived similarly.

We compute a counterpart of Table 2 based on the componentwise 2α/5-quantiles of θ B n, k. They simultaneously provide (1 -2α) = 95%-coverage for parameter M on its own by the applied Bonferroni correction, since M has 5 degrees of freedom, see Table 3.

Elementary algebra permits to compute an evaluation δ(t, ā(t)) of, and confidence lower-bound δ -(t, ā(t)) on, the expected years of life lost for any couple (t, ā(t)) of age t at incident lung cancer and history ā(t) of occupational exposure to asbestos till t. Let us consider three examples:

• Consider a case of incident lung cancer at age t who spent, till that age, 30 years with an occupational exposure to asbestos ε = 332: one evaluates the 95%-confidence lower bound δ -(t, ā(t)) = 30 × 0.004 = 0.09 expected years of life lost (approximately 44 days) and δ(t, ā(t)) = 30 × 0.101 = 3.03 expected years of life lost. This is quite an extreme case, since 3 out of the 8,432 employments described in the dataset achieve the description ε = 332. Table 3: Estimated values (precision 10 -3 ) of the factor of acceleration of time (M (ε) -1) and related right confidence intervals for each level of exposure ε ∈ E \ {0}. Recall that M (0) = 1. A Bonferroni correction ensures that the confidence regions simultaneously guarantee (1-2α) = 95%coverage (for {M (ε) -1 : ε ∈ E} on its own).

ε M (ε) -1 ε M (ε) -1 ε M (ε)
• Consider a case of incident lung cancer at age t who spent, till that age, 10 years (then later 5 years and 2 years) with an occupational exposure to asbestos ε = 322 (then later ε = 121 and ε = 222): one evaluates the 95%-confidence lower bound δ - This is quite an extreme case, since only 6 and 3 out of the 8,432 employments described in the dataset achieve the descriptions ε = 213 and ε = 223.

Among the n = 860 cases of our dataset, only 259 (i.e., 30%) cases are associated with positive expected years of life lost. We report in Table 4 the quartiles, mean and extreme values of expected years of life lost as computed on those 259 cases.

• The maximum value is reached by a male who accumulated through his professional life a total of 33 years with occupational exposure to asbestos equal to ε = 313 and was diagnosed a lung cancer at 70 years old. Although this is not relevant as far as the evaluation of the expected years of life lost is concerned, his lifetime tobacco equals 45 pack years.

• The minimum value is reached by 4 women who accumulated through their professional lives a total of 1 year with occupational exposure to asbestos ε ∈ {211, 221} and were diagnosed a lung cancer at 51 (for two of them), 59 and 68 years old. Although this is not relevant as far as the evaluation of the expected years of life lost is concerned, their lifetime tobacco uses equal 25, 30, 32 and 55 pack years).

• The median value is reached by a man who accumulated through his professional life a total of 4 years (respectively, 5 and 7) with occupational exposure to asbestos equal to ε = 111 (respectively, ε = 211 and ε = 212) and was diagnosed a lung cancer at 71 years old. Although this is not relevant as far as the evaluation of the expected years of life lost is concerned, his lifetime tobacco equals 55 pack years.

We represent in Figure 1 the empirical cdf of the expected years of life lost (and corresponding 95%-confidence lower bounds) for the 259 cases for whom it is positive. min.

25%

50% mean 75% max. expected years of life lost 0.026 0.289 0.769 2.467 2.408 36.577 95%-lower bound 0.001 0.014 0.037 0.555 0.102 12.832 Table 4: Quartiles, mean and extreme values of the expected years of life lost and corresponding 95%-confidence lower-bound (precision 10 -3 ), as computed on those 259 cases (i.e., 30% of all cases) for whom the evaluated expected years of life lost is positive.

years empirical cdf 0.001 0.01 0.1 1 5 10 40 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure 1: Empirical distributions of expected years of life lost and related confidence lower-bound. The rightmost curve with bullets (respectively leftmost curve with triangles) represents the empirical cdf of the expected years of life lost (respectively of the 95%-confidence lower bound on that number) of those cases for whom it is positive, that is the empirical cdf of

{δ(T 1 i , Ā1 (T 1 i )) : δ(T 1 i , Ā1 (T 1 i )) > 0, i ≤ n} (respectively {δ -(T 1 i , Ā1 (T 1 i )) : δ(T 1 i , Ā1 (T 1 i )) > 0, i ≤ n}).
Only 30% of the cases are concerned. The x-axis scale is logarithmic.

Discussion

We have developed a collection of threshold regression models (see Section 4.3), and have dataadaptively selected a better model in it by relying on multi-fold likelihood-based cross-validation (see Sections C.2 and C.3 for the descriptions of the model selection procedure and derived better model). The latter better threshold regression model has been fitted by maximum likelihood, and bootstrapped confidence intervals have been obtained (see Section 5.1). The statistical procedure has been adjusted in order to eliminate the bias induced by the matched case-control sampling design used to collect the dataset (see Sections 3.3 and 4.3). This necessary preliminary step was made possible because the joint distribution of (V, Y ) in the population of interest can be computed beforehand independently from our dataset (see Section C.1). We have discussed the implications of the fitted threshold regression model in terms of expected years of life (free of lung cancer) lost (due to occupational exposure to asbestos) which is naturally attached to it (see Section 5.2).

We finally acknowledge a limitation of the approach undertaken in this article: The link between the occupational exposure to asbestos and age at incident lung cancer is well-defined in the context of the proposed threshold regression models, but we do not extend it beyond. The parameter we aim for is therefore difficult to comprehend (it is related to the Kullback-Leibler projection of the true distribution of the data onto a threshold regression model), and the inference procedure certainly fails to estimate optimally/efficiently what we really care for, which would be a measure of the strength of the link between the occupational exposure to asbestos and age at incident lung cancer defined non-or semiparametrically. We intend to build on the present study and go further in that direction in future work.

A Asymptotic properties of the case-control weighted maximum likelihood estimator

We recall that

F = {ℓ(p ⋆ θ ) : θ ∈ Θ} is P 0 -Glivenko-Cantelli if sup θ∈Θ | 1 n n i=1 ℓ(p ⋆ θ )(O i ) - E P0 ℓ(p ⋆ θ )(O)| = o P0 (1).
The following classical consistency result holds (see Theorem 5.7 and Example 19.8 in [START_REF] Van Der | Asymptotic statistics, volume 3 of Cambridge Series in Statistical and Probabilistic Mathematics[END_REF]).

Proposition 2. Assume that E P ⋆ 0 log p ⋆ 0 (O) is well-defined and that the mapping θ → KL(p ⋆ 0 , p ⋆ θ ) from Θ to the nonnegative real numbers attains its minimum uniquely at

θ 0 ∈ int(Θ) (p ⋆ θ0 is the Kullback-Leibler projection of p ⋆ 0 upon {p ⋆ θ : θ ∈ Θ}). If F is P 0 -Glivenko-Cantelli then θ n converges in probability to θ 0 . This is for instance the case if Θ is a compact metric space, if θ → ℓ(p ⋆ θ )(o ⋆
) is continuous for every o ⋆ and if F admits an integrable envelope function with respect to P 0 .

We also derive an asymptotic normality result (inspired by classical results of asymptotic normality, see Theorem 5.23 in [START_REF] Van Der | Asymptotic statistics, volume 3 of Cambridge Series in Statistical and Probabilistic Mathematics[END_REF]]; we omit the measurability conditions).

Proposition 3 (first part). In the context of Propositions 1 and 2, assume in addition that θ → log p ⋆ θ (Z|Ω) is twice differentiable at θ 0 P ⋆ 0 -almost surely with first and second derivatives l⋆ θ0 (O ⋆ ) and l⋆ θ0 (O ⋆ ) such that l(o ⋆ )dP ⋆ 0 (o ⋆ , Y = y) are properly defined for l = l⋆ θ0 , l⋆ θ0 and y = 0, 1, and introduce accordingly the weighted versions

l(θ 0 )(O) = q 0 l⋆ θ0 (V 1 , O 1⋆ ) + q0 (V 1 ) 1 J J j=1 l⋆ θ0 (V 1 , O 0,j⋆ ), l(θ 0 )(O) = q 0 l⋆ θ0 (V 1 , O 1⋆ ) + q0 (V 1 ) 1 J J j=1 l⋆ θ0 (V 1 , O 0,j⋆ ).
Suppose also that, for every θ 1 , θ 2 in a neighborhood of θ 0 and a function ṁ such that

E P ⋆ 0 ṁ(O ⋆ ) 2 < ∞, P ⋆ 0 -almost surely | log p ⋆ θ1 (Z|Ω) -log p ⋆ θ2 (Z|Ω)| ≤ ṁ(O ⋆ ) θ 1 -θ 2 .
Furthermore, assume that θ → E P0 ℓ(p ⋆ θ )(O) = E P ⋆ 0 log p ⋆ θ (O ⋆ ) (by virtue of Proposition 1) admits a second-order Taylor expansion at θ 0 with nonsingular symmetric second derivative matrix S θ0 =

E P0 l(θ 0 )(O). Then S θ0 = E P ⋆ 0 l⋆ θ0 (O ⋆ ) and √ n(θ n -θ 0 ) = -S -1 θ0 1 √ n n i=1 l(θ 0 )(O i ) + o P0 (1). ( 22
)
In particular, the sequence √ n(θ nθ 0 ) is asymptotically Gaussian with mean zero and covariance

matrix Σ = S -1 θ0 E P0 [ l(θ 0 )(O) l(θ 0 )(O) ⊤ ]S -1 θ0 .
We purposely do not use the same convention to denote the first and second order derivatives at θ 0 of θ → log p ⋆ θ (Z|Ω) (respectively l⋆ θ0 (O ⋆ ) and l⋆ θ0 (O ⋆ )) and the derivatives of θ → l(θ)(O) (respectively l(θ 0 )(O) and l(θ 0 )(O)). We intend to stress that the former are related to the representative sampling observed data structure O ⋆ whereas the latter are related to the casecontrol sampling observed data structure O.

A natural question arises: How does the asymptotic covariance matrix Σ compare with the asymptotic covariance matrix one would have got under representative sampling? We give in the second part of Proposition 3 a very simple answer, but for a representative sampling under a modified version of P ⋆ 0 . Introduce for clarity of exposition the notation q0 (o ⋆ ) = q 0 q 0 (y|v) q 0 (1|v) ( 23)

such that q0 (o ⋆ ) = q 0 if y = 1 (o ⋆ corresponds to a case) and q0 (o ⋆ ) = q0 (v) if y = 0 (o ⋆ corresponds to a control). By setting dP ⋆ 1 dP ⋆ 0 (o ⋆ ) = 1 2q 0 (o ⋆ ) , (24)
we define a probability distribution P ⋆ 1 for the observed data structure O ⋆ (indeed, o ⋆ → q0 (o ⋆ ) is positive and E P ⋆ 0 (2q 0 (O ⋆ )) -1 = 1). Moreover:

• under P ⋆
1 , being a case is as likely as being a control (i.e., P ⋆ 1 (Y = 1) = 1 2 ); • the marginal distribution of the matching variable V under P ⋆ 1 equals the conditional distribution of V under P ⋆ 0 , conditionally on being a case (i.e., P ⋆ 1 (V = v) = q 0 (v|1) for all v ∈ V);

• given (V, Y ), O ⋆ has the same distribution under P ⋆ 1 as under P ⋆ 0 (indeed, q0 (o ⋆ ) depends on o ⋆ through (v, y) only). Furthermore, since obviously

2E P ⋆ 1 q0 (O ⋆ ) log p ⋆ θ (O ⋆ ) = E P ⋆ 0 log p ⋆ θ (O ⋆ ), θ → q0 (O ⋆ ) log p ⋆ θ (O ⋆
) is a proper loss function for the purpose of estimating θ 0 under representative sampling of

O ⋆ ∼ P ⋆ 1 . Proposition 3 (second part). Define S ′ θ0 = E P ⋆ 1 q0 (O ⋆ ) l⋆ θ0 (O ⋆ ).
Suppose that the model is wellspecified (or equivalently that KL(p ⋆ 0 , p ⋆ θ0 ) = 0). Assume in addition that for all o 1⋆ = (z 1 , ω 1 ), the class of derivatives of p ⋆ θ (z 1 |ω 1 ) with respect to θ is uniformly bounded (in θ) by an integrable function (of z 1 ). Then the covariance matrix Σ satisfies

Σ = 1 2 S ′ θ0 -1 E P ⋆ 1 [q 0 (O ⋆ ) 2 l⋆ θ0 (O ⋆ ) l⋆ θ0 (O ⋆ ) ⊤ ]S ′ θ0 -1 .
In particular, 2Σ can be interpreted as the asymptotic covariance matrix of the M -estimator of θ 0 based on the loss function θ → q0 (O ⋆ ) log p ⋆ θ (O ⋆ ) and n iid observations drawn from P ⋆ 1 . The n observations under P 0 -case-control sampling correspond to 2 × n observations under P ⋆ 1representative sampling, each of them in the former counting for two in the latter. Elements of proof are relegated to Section B.

B Elements of proof

Proof of Proposition 1. On one hand, note that

E P0 q 0 log p ⋆ (V 1 , O 1⋆ ) = q 0 log p ⋆ (v 1 , o 1⋆ )dP ⋆ 0 (v 1 , o 1⋆ |y = 1) = log p ⋆ (v 1 , o 1⋆ )dP ⋆ 0 (v 1 , o 1⋆ , y = 1) = log p ⋆ (o ⋆ )dP ⋆ 0 (o ⋆ , y = 1). ( 25
)
On the other hand, for each j ≤ J,

E P0 q0 (V 1 ) log p ⋆ (V 1 , O 0,j⋆ ) = E P0 q0 (V 1 )E P0 log p ⋆ (V 1 , O 0,j⋆ )|V 1 = E P0 q0 (V 1 ) log p ⋆ (V 1 , o ⋆ )dP ⋆ 0 (o ⋆ |V 1 , y = 0) = q0 (v 1 ) log p ⋆ (v 1 , o ⋆ )dP ⋆ 0 (o ⋆ |v 1 , y = 0)dP 0 (v 1 ).
Furthermore, for each v ∈ V, dP 0 (v) = dP ⋆ 0 (v|y = 1) = q 0 (v|1)δ v (v) (we use the same shorthand notation as in ( 4), ( 5), ( 6), ( 7)) and denote by δ v the Dirac mass at v), hence q0 (v)dP 0 (v) = q 0 q 0 (0|v) q 0 (1|v)

q 0 (v|1)δ v = q 0 (0|v)P ⋆ 0 (v)δ v (v) = dP ⋆ 0 (v, y = 0).
Consequently, we obtain

E P0 q0 (V 1 ) log p ⋆ (V 1 , O 0,j⋆ ) = log p ⋆ (v 1 , o ⋆ )dP ⋆ 0 (o ⋆ |v 1 , y = 0)dP ⋆ 0 (v 1 , y = 0) = log p ⋆ (v 1 , o ⋆ )dP ⋆ 0 (v 1 , o ⋆ , y = 0) = log p ⋆ (o ⋆ )dP ⋆ 0 (o ⋆ , y = 0) (26)
(which does not depend on j). Combining ( 25), (26) finally yields

E P0 ℓ(p ⋆ )(O) = log p ⋆ (o ⋆ )dP ⋆ 0 (o ⋆ ) = E P ⋆ 0 log p ⋆ (O ⋆ ).
The conclusion is straightforward, because

E P ⋆ 0 log p ⋆ (O ⋆ ) -E P ⋆ 0 log p ⋆ 0 (O ⋆ ) = -KL(p ⋆ 0 , p ⋆
), the opposite of the Kullback-Leibler divergence between p ⋆ 0 and p ⋆ , which is positive for p ⋆ = p ⋆ 0 and equals zero otherwise.

Proof of Proposition 3. The expansion (22) and the related distributional limit result are a consequence of [van der Vaart, 1998, Theorem 5.23]. The fact that

S θ0 = E P ⋆ 0 l⋆ θ0 (O ⋆ ) is obtained by adapting slightly the proof of Proposition 1. Regarding E P0 [ l(θ 0 )(O) l(θ 0 )(O) ⊤ ], let us abbreviate xx ⊤ to x 2 and note that l(θ 0 )(O) l(θ 0 )(O) ⊤ = q 0 q0 (O 1⋆ ) l⋆ θ0 (V 1 , O 1⋆ ) 2 + q0 (V 1 ) 1 J j q0 (O 0,j⋆ ) l⋆ θ0 (V 1 , O 0,j⋆ ) 2 + q 0 q0 (V 1 ) l⋆ θ0 (V 1 , O 1⋆ ) 1 J j q0 (O 0,j⋆ ) l⋆ θ0 (V 1 , O 0,j⋆ ) ⊤ + q 0 q0 (V 1 ) 1 J j q0 (O 0,j⋆ ) l⋆ θ0 (O 0,j⋆ ) l⋆ θ0 (O 1⋆ ) ⊤ .
The P 0 -expected value of the first term between brackets is E P ⋆ 0 l⋆ θ0 (O ⋆ ), as another simple adaptation of the proof of Proposition 1 straightforwardly yields. Moreover,

E P0 q 0 q0 (V 1 ) l⋆ θ0 (V 1 , O 1⋆ ) 1 J j q0 (O 0,j⋆ ) l⋆ θ0 (V 1 , O 0,j⋆ ) ⊤ = E P0 q 0 q0 (V 1 )E P0 l⋆ θ0 (V 1 , O 1⋆ ) 1 J j q0 (O 0,j⋆ ) l⋆ θ0 (V 1 , O 0,j⋆ ) ⊤ V 1 = E P0 q 0 q0 (V 1 )E P0 l⋆ θ0 (V 1 , O 1⋆ ) V 1 E P0 1 J j q0 (O 0,j⋆ ) l⋆ θ0 (V 1 , O 0,j⋆ ) ⊤ V 1 by conditional independence. Denote by Π = E P0 ( l⋆ θ0 (V 1 , O 1⋆ )|O 1⋆ \ Z 1 ) the conditional expec- tation of l⋆ θ0 (V 1 , O 1⋆ ) given every component of O 1⋆ but Z 1 , that is given Ω 1 (compatible with V 1
). The projection Π can be written as a measurable function of Ω 1 times l⋆ θ0 (z, Ω 1 )p ⋆ θ0 (z|Ω 1 )dz = ∂p ⋆ θ (z|Ω 1 ) ∂θ θ=θ0 dz = 0, provided that the order of differentiation and integration can be reversed. This is ensured by the stated constraint on the derivatives of p ⋆ θ (z|Ω 1 ) with respect to θ. Consequently, the P 0 -expected value of the second term between brackets in the first display is zero, hence the validity of the alternative version of Σ. The conclusion simply follows from another application of [van der Vaart, 1998, Theorem 5.23] in the classical iid framework associated with P ⋆ 1 .

C Application C.1 Conditional distribution of being a case

Estimating the joint distribution of (V, Y ) hence q 0 (4), (q 0 (1|v)) v∈V (5), and (q 0 (v)) v∈V (7) is made possible thanks to [START_REF] Belot | Cancer incidence and mortality in France over the period 1980-2005[END_REF], an independent study of cancer incidence and mortality in France over the period 1980-2005 (for the conditional distribution of Y given V ), and on data made publicly available by the French National Institute of Statistics and Economic studies (for the marginal distribution of V , see http://www.insee.fr/en/). However, we must assume either (i) that the data from [START_REF] Belot | Cancer incidence and mortality in France over the period 1980-2005[END_REF], which are collected over the whole French population, are representative of the Parisian population of interest, or (ii) that sampling from the four Parisian hospitals that participate to the study is stochastically equivalent to sampling from the population of France. We first estimate these quantities for each year from 1999 to 2002 separately. In agreement with our stationarity assumption (1), we remark that the various estimates are very consistent over the years. In order to gain precision, we average the estimates over the years. The final estimates are presented in Table 5. We emphasize that the weights (q 0 (v)) v∈V are far from being homogeneous.

C.2 Multi-fold likelihood-based cross-validation

The likelihood risk of θ ∈ Θ is by definition

R(θ) = -E P0 l(θ)(O),
which is closely related to minus the Kullback-Leibler divergence between the density p ⋆ 0 of P ⋆ 0 and p ⋆ θ , as explained in Section 3. Let us denote by θ n,k (P n ) the case-control weighted maximum likelihood estimator defined in (10) with θ ranging over Θ k . Given the collection {θ n,k (P n ) : k ∈ K} we wish to select the estimator θ n, k(P n ) that minimizes R, where k itself depends on P n . Because q 0 = 470.0682e-06 a q 0 (1|0, a) q 0 (1|1, a) 1 2.058932e-06 1.663324e-06 2 1.859944e-05 1.460473e-05 3 6.803086e-05 4.461827e-05 4 2.586692e-04 1.184914e-04 5 6.484864e-04 1.947058e [START_REF] Belot | Cancer incidence and mortality in France over the period 1980-2005[END_REF]. Left: Estimates of q 0 (1|w 1 , v 2 ), as defined in (5). Middle: Estimate of q 0 , as defined in (4). Right: Estimates of q0 (w 1 , v 2 ), as defined in (7). Here, w 1 = 0 for men and w 1 = 1 for women, and v 2 = a if the age at sampling x belongs to [t a ; t a+1 ), where t 0 = 0, t a = 30 + 5(a -1) for 1 ≤ a ≤ 12 and t 13 = ∞.

the definition of R involves the true distribution P 0 , we must estimate R(θ n,k (P n )) and choose to do so by multi-fold cross-validation. Details follow. We split the data randomly into a training and a validation samples. Given an integer V (later set to V = 10), each observed data structure O i is associated with a label lab i = 1 + (i mod V ).

The collection of labels {lab

i : i ≤ n} ⊂ {1, . . . , V } is such that max l,l ′ ≤V | n i=1 1{lab i = l} - n i=1 1{lab i = l ′ }| ≤ 1.
The splitting random variable S = (S 1 , . . . , S n ) ∈ {0, 1} n is drawn independently of O 1 , . . . , O n in such a way that, for each 1 ≤ l ≤ V , S = (1{lab 1 = l}, . . . , 1{lab n = l}) with probability V -1 . Conditionally on S, the observed data structure O i belongs to the training sample if S i = 0 (there are approximately n(V -1)/V such O i 's), otherwise it belongs to the validation sample. The empirical distribution of those O i 's for which S i = 0 (respectively, S i = 1) is P 0 n,S (respectively, P 1 n,S ). The empirical distribution of those O i 's for which lab i = l (respectively, lab i = l) is P l n (respectively, P -l n ). Each Θ k yields a maximum likelihood estimator θ n,k (P 0 n,S ) based on the training sample only. Its risk, averaged over the splits, writes as crit(k) = E S R(θ n,k (P 0 n,S )) = -

1 V V l=1 E P0 l θ n,k (P -l n )(O) .
The value k that minimizes k → crit(k) over K is called the oracle because it depends both on P n and on P 0 . In our attempt to reach that k which is a good proxy to k, we estimate crit(k) by

C.3 Model selection procedure in action

We explain in Section C.2 how the best model index k (with related best model {P ⋆ θ : θ ∈ Θ k}) is obtained in a pre-determined collection K of sub-model indices (with related sub-models {P ⋆ θ :

θ ∈ Θ k }, k ∈ K). The latter collections are constructed by recursion as presented below.

We first initialize Θ 0 = Θ and K -1 = ∅ with convention max ∅ = 0. At a given step ν ≥ 0, a sub-model Θ ν is defined as a subset of Θ meeting ν independent one-dimensional constraints on M ∈ M (i.e., constraints of the type M k,l-1 = M k,l for some k = 1, 2, 3 and l = 1, 2, 3 with convention M k,0 = 0). Start with c(ν + 1) = -∞ and Θ ν+1 = ∅. The following rule is applied to Θ ν+1 : Rule 1. For every possible Θ ′ ⊂ Θ ν derived from Θ ν by adding another one-dimensional constraint on M as described above (all such models share the same dimension), evaluate the corresponding maximum log-likelihood criterion ℓ(Θ ′ ) = max θ∈Θ ′ P n l(θ). If ℓ(Θ ′ ) ≥ c(ν + 1), update c(ν + 1) = ℓ(Θ ′ ) and Θ ν+1 = Θ ′ .

Applying Rule 1 as long as possible yields 7 sets Θ ν , ν = 0, . . . , 6. Their description is given in Table 6.

Θ 0 = Θ, Θ 1 = {θ ∈ Θ 0 : M 1,1 = 0}, ( 27
) Θ 2 = {θ ∈ Θ 1 : M 2,1 = M 2,2 }, (28)
Θ 3 = {θ ∈ Θ 2 : M 2,2 = 1}, (29)
Θ 4 = {θ ∈ Θ 3 : M 3,1 = M 3,2 }, (30)
Θ 5 = {θ ∈ Θ 4 : M 1,2 = 1}, (31) Θ 6 = {θ ∈ Θ 5 : M 3,1 = 0}.

(32) Table 6: Descriptions of Θ 0 , . . . , Θ 6 . The collection of parameter sets is nested. For instance, Θ 3 is the set of those θ ∈ Θ such that M 1,1 = 0 and M 2,1 = M 2,2 = 1. Regarding dimensions, it trivially holds that, for each 0 ≤ k ≤ 6, dim(Θ k ) = 27k.

At a given step ν ≥ 0, a set K ν-1 of successive integers is defined. Start with K ν = {max K ν-1 + 1} (a set initially containing a single element) and define Θ ν,max K ν-1 +1 = Θ ν . The following second rule is applied to K ν : Rule 2. For every possible constraint "ϕ(θ) = 0" on θ ∈ Θ ν of the form "α and β independent of W l " for some l = 1, 2, 3 (the lth coordinate of W does not affect the value of the initial health and drift parameters h and µ, see ( 19) and ( 20)), update K ν = K ν ∪ {max K ν + 1} and define Θ ν,max K ν = {θ ∈ Θ ν : ϕ(θ) = 0}.

(Note that each Θ ν therefore gives rise to 2 3 = 8 sets Θ ν,l .) We apply Rule 2 for ν = 0, . . . , 6, and finally define K = ∪ 11 ν=0 K ν = {1, 2, 3, . . . , 56}. For every k ∈ K there exists a unique ν = 0, . . . , 7 such that Θ ν,k is defined: setting Θ k = Θ ν,k concludes the definition of the collection {Θ k : k ∈ K} of interest.

The best model {P ⋆ θ : θ ∈ Θ k} (according to our multi-fold likelihood-based cross-validation criterion) is a subset of {P ⋆ θ : θ ∈ Θ 2 }, featuring 16 degrees of freedom. Its complete description follows:

• the initial health parameter depends on W only through gender (hence not on the indicator of occurrence of lung cancer in close family);

• the drift parameter depends on W only through gender and lifetime tobacco use (hence not on the indicator of occurrence of lung cancer in close family);

• exposure to asbestos is significantly noxious; there is no difference between low probability and no exposure to asbestos at all (in view of ( 16), M 1,1 = 0) and no difference either between low and mild frequencies (in terms of ( 16), M 2,1 = M 2,2 ).

Comment on implementation. Our implementation of the model selection procedure relies on an algorithmic trick. The main difficulty arises from the combination of the following two facts:

• in order to fit a given model {P ⋆ θ : θ ∈ Θ ′ } corresponding to a subset Θ ′ ⊂ Θ, it is necessary to provide the optimization algorithm with the set of constraints which characterize Θ ′ ;

• given the huge number of sub-models, it is out of question to prepare beforehand all sets of constraints for all sub-models that we may have to fit in the course of the model selection procedure.

Hence we create an algorithm which maps an explicit description of any Θ ′ ⊂ Θ onto the corresponding set of constraints required by the optimization procedure.

  mean and median values of θ B n, k = {θ n(1-p), k(P b n(1-p) ) : b ≤ B} only very slightly differ from each other. Moreover, they are very close to the maximum likelihood estimate θ n, k(P n ) computed on the whole dataset. The componentwise α/16-and (1α/16)quantiles of θ B n,

  (t, ā(t)) = 10 × 0.004 + 5 × 0 + 2 × 0.004 = 0.048 years of life lost (approximately 17.5 days) and δ(t, ā(t)) = 10 × 0.094 + 5 × 0 + 2 × 0.092 = 1.124 expected years of life lost. Note that 150, 36 and 189 out of the 8,432 employments described in the dataset achieve the descriptions ε = 121, ε = 222 and ε = 322. • Consider a case of incident lung cancer at age t who spent, till that age, 10 years (then later 15 years) with an occupational exposure to asbestos ε = 213 (then later ε = 223): one evaluates the 95%-confidence lower bound δ -(t, ā(t)) = 10 × 0.374 + 15 × 0.374 = 9.350 expected years of life lost and δ(t, ā(t)) = 10 × 1.078 + 15 × 1.078 = 26.95 expected years of life lost.

  and propose to use the value k that minimizes k → crit(k) over K, whose definition is postponed to Section C.3. In conclusion, the final estimator is θ n, k(P n ).

  table of Table1the overall number of employments that feature a particular value of each coordinate of the "probability/frequency/intensity" description.

	ε	nb. of emp.	ε	nb. of emp.	ε	nb. of emp.
	111	213	211	53		311	138
	112	167	212	64		312	105
	113	3	213	6		313	24
	121	150	221	59		321	136
	122	46	222	36		322	189
	123	3	223	3		323	22
	131	0	231	2		331	1
	132	0	232	0		332	3
	133	0	233	0		333	0
				1	2	3
		probability 582 223 618
		frequency 773 644	6
		intensity 752 610 61
	Table

Table 5 :

 5 Estimating the probability distribution of being a case, based on the independent study

		a	q0 (0, a)	q0 (1, a)
		1 228.3063171 282.6071669
		2 25.2727716	32.1855444
		3	6.9091613	10.5348607
		4	1.8167860	3.9666376
	-04	5	0.7243998	2.4137787
	6 1.192778e-03 2.542976e-04	6	0.3936251	1.8480261
	7 1.854668e-03 3.294062e-04	7	0.2529813	1.4265470
	8 2.331553e-03 3.588764e-04	8	0.2011416	1.3093632
	9 2.928415e-03 4.466062e-04	9	0.1600496	1.0520638
	10 3.686216e-03 5.312313e-04	10	0.1270504	0.8843954
	11 3.608302e-03 5.332930e-04	11	0.1298040	0.8809745
	12 3.636995e-03 5.395069e-04	12	0.1287763	0.8708223
	13 2.171286e-03 3.234775e-04	13	0.2160229	1.4527010

  (27) low probability does not differ from no exposure at all; (28) moreover, low and mild frequencies do not differ; (29) moreover, mild and high frequencies do not differ; (30) moreover, low and mild intensities do not differ; (31) moreover, mild and high probabilities do not differ; (32) moreover, low intensity does not differ from no exposure.

We complete the characterization of our model, and derive the related log-likelihood loss function. By Section 3, we know that it suffices to model the distribution of the observed data structure O ⋆ under representative sampling. As explained in Section 3, we wish to model parametrically the conditional distribution of O ⋆ given Ω = (W, X, Ā(X), Y ), i.e., the conditional distribution of Z given Ω, leaving the conditional distribution of Ω given Y unspecified.

* This collaboration took place while A. Chambaz was a visiting scholar at UC Berkeley, supported in part by a Fulbright Research Grant and the CNRS. A. Chambaz would like to thank M-L. Ting Lee for interesting discussions on threshold regression models.

† D. Choudat and J-C. Pairon were partially supported by a grant from the Agence nationale de sécurité sanitaire de l'alimentation, de l'environnement et du travail (ES 2005-006).

‡ The authors would like to thank two anonymous reviewers for their excellent suggestions.