Helena Marti
email: helena.marti-soler@inserm.fr

Michel Chavance
email: michel.chavance@inserm.fr

Multiple imputation analysis of casecohort studies

Keywords:

The usual methods for analyzing casecohort studies rely on sometimes not fully ecient weighted estimators. Multiple imputation might be a good alternative because it uses all the data available and approximates the maximum partial likelihood estimator. This method is based on the generation of several plausible complete data sets, taking into account uncertainty about missing values. When the imputation model is correctly dened, the multiple imputation estimator is asymptotically unbiased and its variance is correctly estimated. We show that a correct imputation model must be estimated from the fully observed data (cases and controls), using the case status among the explanatory variable. To validate the approach, we analyzed casecohort studies rst with completely simulated data and then with casecohort data sampled from two real cohorts. The analyses of simulated data showed that, when the imputation model was correct, the multiple imputation estimator was unbiased and ecient. The observed gain in precision ranged from 8 to 37% for phase-1 variables and from 5 to 19% for the phase 2 variable. When the imputation model was misspecied, the multiple imputation estimator was still more ecient than the weighted estimators but it was also slightly biased. The analyses of case cohort data sampled from complete cohorts showed that even when no strong predictor of the phase-2 variable was available, the multiple imputation was unbiased, as precised as the weighted estimator for the phase2 variable and slightly more precise than the weighted estimators for the phase-1 variables. However the multiple imputation estimator was found to be biased when, because of interaction terms, some coecients of the imputation model had to be estimated from small samples. Multiple imputation is an ecient technique for analyzing case-cohort data. Practically, we suggest building the analysis model using only the case cohort data and weighted estimators. Multiple imputation can eventually be used to reanalyze the data using the selected model in order to improve the precision of the results.

Introduction

Cohort studies, which facilitate causal interpretations, are popular but expensive. Because precision is mainly limited by the number of cases, it is not essential to collect complete information for all the controls. Thus, casecohort studies and nested casecontrol studies enable cost reduction with a minimal loss of eciency [START_REF] Langholz | Nested casecontrol and casecohort methods of sampling from a cohort: a critical comparison[END_REF]. Casecohort studies were initially proposed by Prentice [START_REF] Prentice | A casecohort design for epidemiologic cohort studies and disease prevention trials[END_REF]. When using this approach, the information collected for incompletely observed controls is ignored and inecient estimators for the eect of phase-1 variables are obtained. Also, the weighted estimators used to analyze casecohort data are not fully ecient and this could aect the estimate of the eect of phase-2 variable(s). Alternatively, Breslow et al. [START_REF] Breslow | Using the whole cohort in the analysis of casecohort data[END_REF] suggested optimizing the sampling weights using all the available data. But casecohort studies can also be viewed as a particular example of incomplete data, in which the observation process is controlled by the study organizers. Paik and Tsai [START_REF] Paik | On using the Cox proportional hazards model with missing covariates[END_REF] proposed a simple imputation approach to model censored observations with missing covariates. In the framework of case cohort studies, it would imply the simple imputation of the expected value of the phase-2 variable(s) for incomplete controls. However, that approach ignores the uncertainty concerning imputationmodel parameters and the values to impute according to a given model.

Multiple imputation is a simple and ecient method for analyzing incomplete observations, while taking into account all the levels of uncertainty regarding missing values. For casecohort data analysis, the multiple imputation estimator may provide improved precision, compared to weighted estimators, because it integrally uses the available information and approximates the partial likelihood estimator, which can be more ecient than the weighted estimators.

The objective of this study was to establish multiple imputation as an alternative to weighted analysis of casecohort data. Below, we present the multiple imputation analysis of casecohort studies and validate this approach by comparing its results to those obtained with weighted estimators. First, we used entirely simulated data. Then, we simulated casecohort surveys from two cohorts, the Prospective Epidemiological Study of Myocardial Infarction (PRIME) study, in which no strong surrogate of the chosen phase-2 variable was available, and the National Wilms' Tumor Study (NWTS) data, for which a surrogate was available. For simplicity, we only consider time-constant covariates. Extension to time-varying covariates is considered in the discussion. [START_REF] Prentice | A casecohort design for epidemiologic cohort studies and disease prevention trials[END_REF] Weighted analysis of casecohort studies Casecohort surveys are examples of two-phase designs. First, the cohort is randomly selected from a general population and the phase-1 information is collected for all the subjects. A subcohort is randomly selected and the entire cohort is followed so as to identify the date of occurrence of the event(s) of interest. Then, the phase-2 information, more expensive, is collected for the subcohort subjects and for all the cases, whether or not they belong to the subcohort. Thus, the phase-2 information is not available for controls not belonging to the subcohort. In cohort surveys, where data are available for the whole cohort, the eect of risks factors on the occurrence of events is generally measured by tting a proportional hazards model. In casecohort surveys this model is based on phase-1 and phase-2 variables and the parameters must be estimated from the available incomplete data. In simulations we will consider two phase-1 variables, Z 1 and Z 3 , and one phase-2 variable, Z 2 .

What is lost in terms of eciency, when using a casecohort design rather than a full cohort analysis, can be quantied by the asymptotic relative eciency (ARE). For a casecohort design with simple random sampling it was shown to be [START_REF] Self | Asymptotic distribution theory and eciency results for casecohort studies[END_REF]:

ARE ≈ 1 + 2 1 -α α 1 + 1 -d d log(1 -d) -1 ≈ 1 -γ, ( 1
)
where α is the proportion of the cohort in the subcohort sample, d is the probability of event occurrence and γ is the fraction of missing information. However, when a phase-1 variable is strongly predictive of the phase-2 variable, stratied sampling of the subcohort can improve eciency as compared to simple random selection [START_REF] Borgan | Exposure stratied casecohort designs[END_REF].

Weighted estimators of the log-relative risks maximize a weighted pseudolikelihood ( L(β)):

L(β) = j      exp{β Z i }w i j k C∪D Y k (t j )exp{β Z k }w k j      (2)
where event j occurs at time t j , C is the subcohort of size n sc , D is the set of cases, Y k (t j ) indicates whether subject k is at risk at time t j , β is the vector of log relative risks, Z k the vector of covariates for subject k, w k j the weight of subject k at time t j , i j the index of the subject whose event occurs at t j , and the symbol denotes transposition. Barlow [7] proposed weighting each complete observation by the inverse of its probability of being included (1 for the cases). Other authors have proposed variable weights, as a function of time, to slightly improve eciency [START_REF] Borgan | Exposure stratied casecohort designs[END_REF].

The variance of this estimator must take into account the increased uncertainty associated with the randomized selection of the subcohort. This requirement can be achieved by using the sandwich variance, which can be estimated as:

V ar( β) = I -1 + n sc (n -n sc ) n CovD C ( 3
)
where I is the Fisher information matrix, n and n sc are the respective sizes of the cohort and subcohort, and CovD C is the empirical covariance matrix of dfbeta residuals from subcohort members dened as [START_REF] Langholz | Computational methods for casecohort studies[END_REF]:

df beta ji = β j -β j(i) s (i) (Z Z) -1 jj (4)
with β j(i) the parameter j estimate obtained after deletion of subject i and s (i) the standard error of this estimate. Stratied sampling of the subcohort considers some information obtained during phase-1, but the information provided by the phase-1 variables is generally ignored in the analysis. Kulich and Lin [START_REF] Kulich | Improving de eciency of relative-risk estimation in case cohort studies[END_REF] proposed a family of doubly weighted estimators intended to more eciently account for the information provided by the initial variables. Qi et al. [START_REF] Qi | Weighted estimators for proportional hazards regression with missing covariates[END_REF] developed nonparametric methods to estimate selection probabilities and nonparametric kernel-smoothing techniques to estimate conditional expectation in fully augmented weighted estimating functions. Breslow et al. [START_REF] Breslow | Using the whole cohort in the analysis of casecohort data[END_REF] suggested calibrating or estimating the weights using all the phase-1 information in order to improve precision: 1) with calibration, the weights are subjected to the constraint that the cohort totals of some auxiliary variables are equal to their weighted sum among all phase-2 subjects. Practically, one builds a prediction model for the phase-2 variable to perform a simple imputation of the predicted values among the controls not belonging to the subcohort, ts the model of interest to the completed data set and uses the inuence function from the model of interest to calibrate. Eventually, the model of interest is tted to the calibrated case-cohort data; 2) with estimation, the weights are the reciprocals of the inclusion probabilities, as estimated from a logistic model tted to the full cohort.

Incomplete observations and multiple imputation

Little and Rubin [START_REF] Little | Statistical Analysis with Missing Data[END_REF] distinguished three observation processes: data missing completely at random (MCAR), when the probability of incomplete observation is constant; data missing at random (MAR), when this probability depends only on observed values; and data missing not at random (MNAR), when this probability depends on unobserved values. The distinction between MAR and MNAR is of utmost importance, because with the former, unbiased estimators of the parameters of interest are available. Casecohort data are MAR, because the probability of being completely observed depends only on case status and, under stratied sampling, on some phase-1 variables.

The multiple imputation method developed by Little and Rubin [START_REF] Little | Statistical Analysis with Missing Data[END_REF] provides an approximation of the maximum likelihood estimator and thus enables the potential selection bias to be corrected. This method relies on the generation of several plausibly completed data sets (M ≥ 2), accounting for all the levels of uncertainty concerning the missing values. A prediction model must be built, taking into consideration the relationships between the incomplete variable and the other variables, as observed in the complete part of the data. The missing data are not replaced by their expectation but by a value drawn from the distribution posited by the model. To take into account the uncertainty concerning the parameters of the imputation model, several imputations are performed with parameters drawn from the asymptotic distribution of their estimator. An estimate of the parameter of interest, θm , m = {1, . . . , M }, and an estimate of the variance of the estimator, V ( θm ), are obtained from each completed data set. If the imputation model is correct, these estimators are not biased. The multiple imputation estimate, also unbiased, is the mean of these M estimates:

θIM = 1 M M m=1 θm (5)
The multiplicity of imputations enables a correct estimation of the variance of this single estimator. The variance is the sum of two components: the withinimputations component (W IM ), estimated as the mean of the M asymptotic variances, W IM , and the between-imputations component (B IM ), estimated from the observed variance of the M estimates, B IM :

V ( θIM ) = W IM + B IM = 1 M M m=1 V ( θm ) + (1 + M -1 ) M m=1 ( θm -θIM )( θm -θIM ) M -1 (6)
where the factor 1 + M -1 is an adjustment for using a nite number of imputations [START_REF] Rubin | Multiple imputation for interval estimation from simple random samples with ignorable nonresponse[END_REF]. Rubin [START_REF] Rubin | Multiple Imputation for Nonresponse in Surveys[END_REF] showed that the relative eciency of multiple imputations with a nite number M , as compared to an innite number of imputations, is:

ARE ≈ 1 + γ M ( 7
)
where γ is the fraction of missing information:

γ ≈ B IM B IM + W IM (8)
A nonexhaustive review of published casecohort studies showed that, among 25 studies, the fraction of missing information ranged from 0.05 to 0.5, with a median around 0.3. With as much as 40% information missing, M=5 imputations provides an ARE=0.97, and, with 50% missing information, M=10 provides an ARE=0.98. Thus, in most instances, there is not much to gain by using more than 5 or 10 imputations to analyze casecohort data.

Multiple imputation requests a correct model of the relationships between the incomplete variable(s) and the variables that are linked to the former. When the statistician doing the multiple imputation is independent of the statistician conducting the analyses: ". . . it is important to include as predictors as many of the variables that are likely to be used in subsequent analyses as possible. Leaving out such variables, even when they are weak predictors, implies that it is known with certainty that they have no relation with the missing values. The result is that correct uncertainty is not reected [START_REF] Rubin | Multiple imputation for interval estimation from simple random samples with ignorable nonresponse[END_REF]." In casecohort studies, the imputation model and the analysis are in the hands of the same statistician, thus only variables useful for the analysis of interest have to be included.

We need to impute missing phase-2 variable values for the controls who do not belong to the subcohort. This requires an imputation model taking into account the dierences between cases and controls: otherwise, the multiple imputation estimator of the eect of the phase-2 variable would be biased. Under the rare disease assumption, it can be shown that a simple generalized linear model using all the casecohort data and including the status indicator among the explanatory variables has to be considered.

Let us assume that the distribution of Z 2 belongs to the exponential family and depends on a phase-1 variable, possibly multidimensional z2 through:

f (z 2 | z2 ) = exp θz 2 -b(θ) + c(z 2 ) a(φ) , ( 9
)
where φ is the dispersion parameter and the canonical parameter θ is a linear function of unknown parameters:

θ = α 0 + α 1 z2 . ( 10
)
Under the rare disease assumption, the distribution of Z 2 is approximately the same for the whole population and among controls:

f (z 2 | z2 , ∆ = 0) f (z 2 | z2 ), ( 11
)
where ∆ is the case indicator. Let π(z 2 , µ 2 , t c ) be the probability of being a case at the end of the observation time, for a subject with Z2 = z2 ,

Z 2 = µ 2 = E[Z 2 | Z2 = z2 ] and censoring time T C = t c ; let π(z 2 , z 2 , t c
) be the probability of being a case at the end of the observation time, for a subject with Z2 = z2 , Z 2 = z 2 and T C = t c . According to the proportional hazards model and using the rare disease assumption:

π(z 2 , z 2 , t c ) = π(z 2 , µ 2 , t c )exp[β(z 2 -µ 2 )] (12)
We also assume that the distribution of the censoring time is independent of Z 2 , so integrating over the censoring time:

π(z 2 , z 2 ) = π(z 2 , µ 2 )exp[β(z 2 -µ 2 )] (13)
and the distribution of Z 2 conditionally on being a case and on Z2 , can be obtained as:

f (z 2 | ∆ = 1, z2 ) = P [∆ = 1 | z2 , z 2 ] P [∆ = 1 | z2 ] f (z 2 | z2 ) = π(z 2 , µ 2 )exp[β(z 2 -µ 2 )] π(z 2 , µ 2 )exp[β(z -µ 2 )]f (z | z2 )dz f (z 2 | z2 ) = exp(βz 2 ) exp(βz)f (z | z2 )dz f (z 2 | z2 ) (14)
The denominator can be developed as:

exp(βz)f (z | z2 )dz = exp βz + θz 2 -b(θ) + c(z 2 ) a(φ) dz = exp (θ + a(φ)β)z -b(θ + a(φ)β) + b(θ + a(φ)β) -b(θ) + c(z) a(φ) dz = exp b(θ + a(φ)β) -b(θ) a(φ) (15)
with θ given by [START_REF] Qi | Weighted estimators for proportional hazards regression with missing covariates[END_REF]. The results of ( 15) and ( 11) lead to:

f (z 2 | ∆ = 1, z2 ) exp a(φ)βz 2 -b(θ + a(φ)β) + b(θ) a(φ) f (z 2 | ∆ = 0, z2 ) = exp [θ + a(φ)β]z 2 -b[θ + a(φ)β] + c(z 2 ) a(φ) (16)
Thus, under the rare disease assumption, the distributions of Z 2 , given ∆ = 0 or = 1, dier by a shift of a(φ)β on the scale of the canonical link, i.e., a shift of σ 2 β for a linear model, a shift of β on the logit scale for a logistic model, or a shift of β on the log scale for a log-linear model, where β is the coecient associated with Z 2 in the proportional hazards model of interest.

The building of the prediction model, particularly the choice of the variable(s) z2 , is crucial to perform multiple imputation. Practically, in addition to the status indicator and stratication variables, it is necessary to adjust for the confounding variables included in the Cox model and for other predictive variables, which could be available.

The analyses were performed with R software (version 2.9.0, The R Foundation for Statistical Computing), using the mice (Multivariate Imputation by Chained Equations) package http://web.inter.nl.net/users/S.van.Buuren/mi/hmtl/mice.htm), which generates multiple imputations, and the survival package, which enables casecohort designs to be carried out and analyzed by means of weighted estimators.

Validation of the method

To validate the method, rst, we used entirely simulated data and compared the estimates and their standard errors to the true values. Then, using the PRIME cohort data, for which no strong surrogate for the chosen phase-2 variable was available, and the NWTS cohort data, which had an available surrogate, we compared the multiple imputation estimator to three weighted estimators: inverse probability weights (the most popular one), calibrated weights and re-estimated weights.

Simulations

Completely simulated data

Two phase-1 variables were simulated: a binary variable, Z 1 , and a Gaussian variable, Z 3 , observed in the entire cohort. Also simulated was a phase-2 standard Gaussian variable, Z 2 , which was independent of Z 1 , but had a correlation of 0.2 with Z 3 . The time to the event of interest had an exponential distribution, with

λ = exp(β 1 Z 1 + β 2 Z 2 + β 3 Z 3 ). β 1 ,
β 2 and β 3 were xed at the same value of 0 or log(2). The censoring time followed a uniform distribution over the interval [0, τ ], where τ was dened so that the probability of the event was approximately 0.01 (τ = 0.008). The cohort size was 25,000 subjects. We also simulated a phase-1 variable predictive of the variable Z 2 , Z 2 = Z 2 + ε with ε ∼ N (0, 1) independent of Z 2 (the correlation between Z 2 and Z 2 was √ 2/2 0.7). We wanted to estimate the eect of Z 2 on the occurrence of the event, adjusting for Z 3 within the framework of stratied sampling of the subcohort. The cohort was divided into 9 strata based on Z 2 and Z 3 tertiles, and the controls were chosen by stratied sampling. Casecohort sampling was simulated with 1,000 subjects in each subcohort (Table 1).

We built a linear prediction model for Z 2 based on the stratum indicators (phase-1 variables) and the status indicator. Z 3 was not directly included in the imputation model to predict Z 2 , because it was used to dene the strata, included in the model, and weakly correlated with Z 2 .

The imputation model was:

Z 2 = α 0 + α 1 status + α 2 I strata + e
where I strata is the vector of stratum indicators. The mean multiple R 2 was: 0.40.

Five plausible sets of complete data were generated for each cohort. Using weighted analysis, the variance estimator II proposed by Borgan et al. [START_REF] Borgan | Exposure stratied casecohort designs[END_REF] was used. One thousand cohorts were simulated.

To assess the consequences of a misspecication of the imputation model, we compared the multiple imputation estimates, obtained when a predictive variable was omitted or included, to weighted estimates. Two scenarios were considered. In the rst, the omitted variable was a confounder. Data were simulated according to the previously used conditions. The cohort was stratied only according to Z 2 tertiles, ignoring Z 3 . We used two imputation models, including the status and strata indicators with or without Z 3 . The same two models were used in the simple imputation stage of the calibrated weighted analysis. For the second scenario, the omitted variable was not a confounder. Data were still simulated according to the previous conditions. The subcohort was selected by simple random sampling so that the imputation model included no strata indicator; Z 2 , which was related to the phase-2 variable (correlation ρ = 0.7), was omitted or included in the model.

Results

In Table 2, for phase-1 (Z 1 and Z 3 ) and phase-2 (Z 2 ) variables, all estimates of the log relative risks were unbiased. Likewise, the multiple imputation variance and Borgan's variance estimator II (BII) agreed with the observed dispersion of estimates. For phase-1 variables, this observed dispersion was close to those of the entire cohort and multiple imputation but larger with the weighted estimator. For the estimated eect of the phase-2 variable, the observed standard errors were obviously larger for casecohort than full cohort analyses, but they were slightly smaller with multiple imputation than with weighted estimators. For phase-1 variables, the relative increases of standard errors for the weighted analysis estimators compared to those of multiple imputation, ranged from 8 to 37%, whereas they were slightly smaller (5 or 19%) for phase-2 variables.

Table 3, gives the results obtained with the calibrated weighted estimator and the multiple imputation estimator using a correct imputation model and a misspecied imputation model omitting a confounder variable. All calibrated weighted estimates were unbiased, for phase-1 and phase-2 variables, using a correct or a misspecied simple imputation model. All multiple imputation estimates were unbiased when a correct imputation model was used. However, the misspecied model yielded biased multiple imputation estimates for the phase-2 variable and the omitted phase-1 variable. All multiple imputation estimates were more precise than calibrated weighted estimates, for phase-1 and phase-2 variables eects. The relative increase of the standard error of the calibrated weighted estimate, as compared to multiple imputation estimate, exceeded 10%.

Table 4 reports the results obtained with the weighted estimator, a correct imputation model and a misspecied imputation model omitting a predictive variable that was not a confounder. All the estimates of phase-1 variable eects were unbiased. The precision of the two multiple imputation estimates, obtained with the correct or the misspecied models, were similar. The weighted estimate was less precise than multiple imputation estimates, more specically, the relative increase of the standard error for the weighted estimator compared to multiple imputation estimators was 20% for Z 1 and over 46% for Z 3 .

For the phase-2 variable, the multiple imputation estimate with the correct imputation model was unbiased, while the misspecied imputation model and the weighted analysis led to estimates slightly biased in opposite directions, respectively, 2.3% and 1.9%. As expected, ignoring Z 3 in the imputation model decreased the precision of the multiple imputation estimator. The relative increase of the standard error with the misspecied model compared to that with the correct model was 7.6%. The standard error of the weighted estimator was larger than that of the multiple imputation estimator. The relative increase of the standard error of the weighted estimator, compared to multiple imputation estimator, was 30% with the correct model and 21% with the misspecied model.

Description of the data

The PRIME survey [START_REF] Yarnell | The PRIME study: classical risk factors do not explain the severalfold dierences in risk of coronary heart disease between France and Northern Ireland[END_REF] was a multicenter cohort investigation studying risk factors of ischemic heart disease (IHD) and other cardiovascular end points. Among the 9,520 male subjects, 642 (6.7%) experienced a cardiovascular event. The median follow-up time was 10 years and for those who suered an event, the median time to its occurrence was 5 years.

We chose to estimate the eect of brinogen (phase-2 variable) on the occurrence of the event, adjusting for age, center, total cholesterol (CH), high-density lipoprotein cholesterol (HDL), systolic blood pressure (SBP) and tobacco use. Casecohort data were simulated based on the entire PRIME cohort, and the results obtained with the dierent estimators were compared to each other and to the results obtained from the full cohort. The validation procedure used 1,000 simulated subcohorts of size 2,100.

Sampling of the subcohort

Fibrinogen is a protein whose circulating concentration increases during inammatory conditions. It also plays an important role in normal and pathological blood coagulation, and elevated brinogen levels are associated with a higher risk of cardiovascular events. Because smoking is one of the main factors determining the brinogen level [START_REF] Scarabin | Plasma brinogen explains much of the dierence in risk of coronary heart disease between France and Northern Ireland. The PRIME study[END_REF], the cohort was stratied according to tobacco use, treated as a phase-1 variable available for the entire cohort. The strata were created as follows: stratum 1: non-smokers, stratum 2: ex-smokers, stratum 3: smokers [19] g/day, stratum 4: smokers [1019] g/day, stratum 5: smokers ≥ 20 g/day. Stratied sampling was used to select subcohorts and the number of events in each stratum is given in Table 5. The sampling probabilities were approximately 0.16, 0.20, 0.21, 0.25 and 0.46, respectively, for strata 15.

The linear model for brinogen imputation used as explanatory variables the status indicator and the variables included in the Cox model: tobacco, age, center, CH, HDL and SBP (mean multiple R 2 : 0.07). Five imputations were performed for each subcohort.

Results

The multiple imputation estimator was compared to the standard weighted estimator [START_REF] Borgan | Exposure stratied casecohort designs[END_REF], calibrated weights estimator [START_REF] Breslow | Using the whole cohort in the analysis of casecohort data[END_REF] and re-estimated weights estimator [START_REF] Breslow | Using the whole cohort in the analysis of casecohort data[END_REF]. The mean of the 1,000 log relative risk estimates, corresponding to the 1,000 subcohorts, are given in Table 6. The median fraction of missing information about the brinogen eect was 0.22, so ve imputations can be considered sucient. For the phase-2 variable, the estimates were similar but, as expected, the standard error was larger in the casecohort analysis than the full cohort analysis. Similar mean standard errors were obtained with multiple imputation and standard weighted estimation (Borgan's variance estimator II), they were very close to those obtained with calibrated and re-estimated weights. For phase-1 variables, the multiple imputation estimates were more precise than those obtained by weighted, calibrated and re-estimated weights analyses, and were nearly the same as those obtained from the entire cohort.

NWTS data

The NWTS cohort [START_REF] Angio | Treatment of Wilms' tumor. Results of the Third National Wilms' Tumor Study[END_REF] consisted of 3,915 patients with Wilms' tumor diagnosed during 19891994 and followed until the earliest sign of disease progression or death for event-free survival. Baseline covariates included stage (IIV), age at diagnosis, tumor diameter and two binary histological evaluations (favorable vs. unfavorable): the local hospital histology and central histology evaluated in a centralized reference laboratory. The former was strongly predictive of the latter (specicity 98%, sensitivity 74%) and both were available for all the patients. However, like Breslow et al. [START_REF] Breslow | Using the whole cohort in the analysis of casecohort data[END_REF], we simulated casecohort studies using central histology as the phase-2 variable. For the NWTS analysis, the Cox model included central histology (phase-2 variable), age as a piecewise linear variable with change point at 1 year, stage III/IV versus I/II, tumor diameter and the interactions local histology*age and stage*diameter. Breslow et al. [START_REF] Breslow | Using the whole cohort in the analysis of casecohort data[END_REF] dened 16 strata according to event-free survival, stage (I/II or III/IV), favorable local histology (Yes or No) and age < 1 year (Yes or No). They sampled controls from only three strata, dened by favorable local histology, and "age < 1 year + stage I/II" (n = 120), "age ≥ 1 + stage I/II" (n = 160) and "age ≥ 1 + stage III/IV" (n = 120), while including all the subjects in the 13 other strata. They predicted unfavorable central histology, according to local histology, stage, age, tumor diameter and the interaction local histology*stage.

These data present two specic features: rst, a phase-1 variable, strongly predictive of the phase-2 variable, is available and, second, an interaction exists between the phase-2 variable and a phase-1 variable. Among these real data, in the sampled strata (all with favorable local histology) only a few controls had unfavorable central histology, especially in stratum 1 (n = 2) and stratum 4 (n = 17) (Table 7). Thus, for some subcohorts, the specic imputation model for these strata could not be estimated because of observed innite odds ratio. Moreover, even when the subcohorts included some controls with unfavorable central histology in these strata, their number was necessarily small, and the estimator of the imputation model parameters might not have been distributed as assumed according to the asymptotic results.

Because of both particularities of these data, we considered two imputation models. The rst, based on the proportional hazards model of interest and local histology, included the status indicator, local histology, stage, age, tumor diameter and the interaction local histology*stage; the second was limited to the status indicator and local histology. When the imputation model included the interaction local histology*age, biased estimates and large standard errors were observed due to the small number of unfavorable central histologies in some subgroups used for the estimation of the interaction terms (Table 8). For the imputation model including only the status indicator and the surrogate variable, the estimates were slightly biased. The simple eect of central histology, which represented the expected difference at age 1 between children with favorable or unfavorable central histology, diered slightly from the eect estimated for the full cohort (respectively, 4.11 vs 4.04); the same held true for the age eect after age 1 for the children with favorable central histology (respectively, 0.10 vs 0.11), and for the age eect after age 1 for the children with unfavorable histology (respectively, 3.63 vs 3.30). With the second imputation model, including only local histology and status indicator, we observed a small bias and good precision: although the multiple imputation was slightly biased, in particular for the age*histology interaction, it was always more precise than the standard weighted estimator and slightly more precise than the calibrated and re-estimated weighted estimators.

Discussion

The aim of the weighted analysis and of multiple imputation is to reconstitute the whole cohort. The former, weights the subjects in the casecohort sample by the inverse of the probability of being observed during phase-2, but the phase-1 data observed for the other subjects are generally ignored. Alternatively, Breslow et al. (2009), proposed two approaches using all the phase-1 information. Multiple imputation uses all the available data supplementing them with plausible values agreeing with what is observed in the complete data.

A key aspect of multiple imputation is the construction of the prediction model.

It is necessary to reproduce correctly the relationship between the outcome and the incomplete variable, adjusting for the confounders included in the Cox model. With casecohort data, the problem is complicated by the censoring process. One might think that it useful to include censoring time in the imputation model because, when the phase-2 variable is predictive of the event, its distribution among controls might not be the same at the beginning and the end of follow-up. However, in section 3, we demonstrated that, for a phase-2 variable with an exponential family distribution, a generalized linear model, including the case-status indicator as an exploratory variable, provides an approximately unbiased multiple imputation estimate. The proof relies on several assumptions: rst, the studied event has to be rare; and second, the phase-2 variable has to be independent of the censoring time. The rst assumption is precisely what justies the use of a case cohort design, while the second is required by the proportional hazards model. Thus, they do not represent new limitations. Our simulations showed that using all the complete subjects, cases as well as controls, and with a correct imputation model, including the case-status indicator, the multiple imputation estimator was unbiased. We also performed some simulations that conrmed no improvement of the multiple imputation estimator by adding the censoring time to the imputation model (data not shown). Confounding variables appearing in the analysis model also have to be included in the imputation model, if they remain predictive of the phase-2 variable, adjusting for the other predictors. On the other hand, misspecication of the imputation model can aect the phase-2 variable estimates but also the estimate of the omitted phase-1 variable. Omitting variables that improve the prediction can yield biased estimates and increase uncertainty about the parameter. Including variables that do not improve the prediction increases the uncertainty of the model coecients and thus the between-imputation variance. Calibrated weighted estimates were found to be less sensitive to a misspecication of the imputation model than multiple imputation estimates concerning bias. However, multiple imputation estimates were more precise. The completely simulated data showed that, when the imputation model was correct, the multiple imputation approach provided unbiased and ecient estimators. Both the weighted and multiple imputation estimators were centered on the true values. An important result of this simulation study was that multiple imputation correctly estimated the variance of its estimators (just as BII correctly estimated the variance of the standard weighted estimator). This nding allowed us to compare the variance estimators of casecohort data simulated from cohort surveys, for which the variance of the estimators cannot be compared to the true value because it is unknown. As expected, the multiple imputation approach was more precise than the usual weighted estimators for the parameters associated with phase-1 variables. The former also was slightly more precise than the latter for the phase-2 variable, despite the fact that it only used a categorized transformation of the explanatory variables (the stratum indicators used for stratied sampling). One explanation might be that multiple imputation approximates the maximum partial likelihood estimate, which is more ecient than weighted estimators. The simulations implying misspecied imputation models revealed, as expected, that the omission of a confounding variable from the imputation model had consequences in terms of bias and precision. We insist that is essential to include the stratication variable(s) and the variables included in the analysis model in the imputation model. The consequences of omitting variable(s) related to the phase-2 variable, but not to the event of interest, mainly concern precision. In these simulations, the weighted estimators did not suer from serious bias problems, not even the calibrated weighted estimator using a misspecied model in its imputation phase. However, they suered from appreciable losses of precision. It should be underlined that most casecohort surveys are generally performed to answer several scientic goals dealing with dierent diseases and exposures. Thus, it can be dicult to dene a stratied sampling ecient for all the analyses and to optimize the weighted estimators. By contrast, multiple imputation can be adapted to each phase-2 variable and each sub-study to improve the precision of the estimates of interest. When a misspecied imputation model was used, multiple imputation estimators for the phase-1 variable eects were not biased and their precision was similar to that obtained with a correct imputation model. For the phase-2 variable eect, slightly biased estimates were observed with the misspecied imputation model and weighted analyses. The eect of the misspecication on the imputation model was more noticeable in terms of precision. The loss of precision using a misspecied model as compared to a correct model was 7.6%. This loss was greater for the weighted estimate than the multiple imputation estimate, and exceeded 21%. We did not include the Z 3 variable, which was weakly correlated to the phase-2 variable. Results were less satisfactory when the correlation between the two variables was stronger, in particular concerning the precision of the weighted estimators (data not shown).

The casecohort data simulated from the PRIME cohort study showed that multiple imputation can be used, even when no strong predictor of the phase-2 variable is available. Using the variables of the analysis model plus the case indicator in the imputation model, the multiple imputation estimator was more precise than the weighted estimators, particularly the standard estimator, for the eects of phase-1 variables. For the phase-2 variable, the multiple imputation estimator had the same precision as the standard weighted estimators and it was slightly less precise than calibrated and re-estimated weighted estimators.

The NWTS cohort data represent one of the worst possible situations for using multiple imputation. Inclusion of an interaction term between the phase-2 variable and the two age-eect components required an imputation model including similar interactions between the indicator status and the age eects. The corresponding coecients had to be estimated in separate strata, with very low numbers of patients presenting unfavorable central histology or even no such patient at all. As a consequence, the maximum likelihood estimator of the imputation model could be expected to be biased [17, chap. 8.4] and the multiple imputation estimator reected this bias. Although local histology was strongly predictive of central histology, imputation model 2, using only the former and the case indicator, also yielded biased estimates of the Cox model parameters. The imputed values correctly reected the global proportion of patients with unfavorable central histology, but not always the proportional dierence between cases and controls, as was the case for the interaction terms age*local histology or stage*diameter.

It is reasonable to wonder how many imputations are needed. The number of requested imputations increases with the proportion of missing information. However, in casecohort studies, the proportion of missing information is considerably smaller than the percentage of incompletely observed subjects, and a small number of imputations, 510, should suce. With the PRIME data and using a small sampling rate (700/25,000), the results obtained with 5 imputations did not dier appreciably from those obtained with 10 (data not shown).

Herein, we presented simulations with only one phase-2 variable. However, the approach can easily be extended to several phase-2 variables. If several covariates are incomplete, we suggest imputing them using a multivariate distribution, which takes into account the correlation structure between the covariates, for instance with the mice package, which generates multivariate imputations via Markov Chain Monte Carlo (MCMC) algorithm according to their joint distribution.

We focused on situations in which covariates were time-xed. However, the multiple imputation approach can be extended to time-varying covariates by, using mixed models to account for the repeated measurement, and the within-subject correlation structure [START_REF] Yucel | Multiple imputation inference for multivariate multilevel continuous data with ignorable non-response[END_REF].

When the phase-2 data allow consistent estimation of the imputation model, multiple imputation is an ecient technique to analyze the data from casecohort surveys. For phase-1 variables, multiple imputation has better eciency than weighted estimators, a valuable improvement for prediction studies or when the eect of some phase-1 variables is a focus of interest. A large number of imputations is not required to obtain good quality estimates. Software that simply implements this procedure is available: under R, the mice library; under Stata, the Imputation by Chained Equations library, or under SAS, the PROC MI and PROC MIANALYZE.

To gain time and determine whether multiple imputation provides estimates similar to weighted analysis, we suggest building the analysis model using only the casecohort data and weighted estimators. Multiple imputation can eventually be used to reanalyze the data with the selected nal model to improve the precision of the results.

: Z2i = α0 + α1Indcas + α2Strata + α3Z3 + ei (correct model) IM2: Z2i = α0 + α1Indcas + α2Strata + ei (misspecied model) Calibrated1: Z2i = α0 + α1Indcas + α2Strata + α3Z3 + ei (correct model) Calibrated2: Z2i = α0 + α1Indcas + α2Strata + ei (misspecied model).
: Z2i = α0 + α1Indcas_i + α2 Z2 + α3Z3 + ei (correct model) IM2: Z2i = α0 + α1Indcas_i + α3Z3 + ei (misspecied model).

Table 1 :

 1 Distribution of subcohort by stratum.

	Correlation (Z2, Z2)=0.7

Table 2 :

 2 Parameter estimates, stratied sampling of the subcohort (mean results from 1,000 simulations). Z 2i = α 0 + α 1 Ind cas_i + α 2 Strata + e i .

	Correlation (Z2, Z2)=0.7

Table 3 :

 3 Consequences of a misspecication of the imputation model. Stratied sampling of the subcohort. Results from 1,000 simulations.

		Cohort	IM1	IM2	Calibrated1	Calibrated2
		β	β	β	β	β
	Parameter	( SE)	( SE)	( SE)	( SE)	( SE)
	β1 = 0.6931	0,7133	0,7017	0,6947	0,7126	0,7123
		(0,1737)	(0,1742)	(0,1743)	(0,1901)	(0,1933)
	β2 = 0.6931	0,6940	0,6867	0,7554	0,6921	0,6902
		(0,0588)	(0,0718)	(0,0694)	(0,0853)	(0,0865)
	β3 = 0.6931	0,6955	0,6911	0,7543	0,6975	0,6970
		(0,0576)	(0,0596)	(0,0573)	(0,0702)	(0,0772)
	Mean of the 1,000 estimates (mean of the 1,000 standard error estimates)
	IM1

Table 4 :

 4 Simple random sampling of the subcohort. Parameter estimates with correct and misspecied imputation model. Mean results from 1,000 simulations SE, mean of the standard error estimates; SE, standard error of the estimates; IM1

	Correlation (Z2, Z2)=0.7

Table 5 :

 5 Distribution of cases by stratum.

	Stratum	Cases (%) Stratum size Subcohort size
	1 Non-smokers	153 (5.3)	2,890	475
	2 Ex-smokers	261 (6.4)	4,078	800
	3 Smokers 1-9 g/day	51 (6.3)	816	175
	4 Smokers 10-19 g/day	55 (7.8)	707	175
	5 Smokers ≥ 20 g/day 122 (11.9)	1,029	475
	Total	642 (6.7)	9,520	2,100

Table 6 :

 6 Estimates of the log relative risks. Unique estimates provided by the PRIME cohort.2 Mean estimations of the 1,000 subcohorts.3 Asymptotic standard error (SE) of the estimate. FN, brinogen; CH, cholesterol; HDL, high-density lipoprotein; SBP, systolic blood pressure. Imputation model: F N = β0 + β1status + β2tobacco + β3age + β4center + β5CH + β6HDL + β7SBP +

		Whole cohort 1	Multiple imputation 2	Standard 2	Weighted estimator Calibrated 2 Re-estimated 2
	Variable	β (SE)	β (SE 3 )	β (SE 3 )	β (SE 3 )	β (SE 3 )
	Fibrinogen	0.1312	0.1381	0.1346	0.1338	0.1337
		(0.0346)	(0.0421)	(0.0425)	(0.0414)	(0.0413)
	Tobacco	0.0219	0.0218	0.0220	0.0214	0.0217
		(0.0036)	(0.0037)	(0.0045)	(0.0040)	(0.0039)
	Age	0.0573	0.0566	0.0573	0.0566	0.0570
		(0.0138)	(0.0138)	(0.0160)	(0.0162)	(0.0160)
	Center	0.2788	0.2764	0.2795	0.2749	0.2770
		(0.0857)	(0.0863)	(0.1006)	(0.1001)	(0.0990)
	CH	0.0078	0.0077	0.0078	0.0078	0.0078
		(0.0010)	(0.0010)	(0.0012)	(0.0012)	(0.0012)
	HDL	-0.0526	-0.0528	-0.0528	-0.0529	-0.0529
		(0.0070)	(0.0070)	(0.0080)	(0.0084)	(0.0082)
	SBP	0.0104	0.0104	0.0105	0.0105	0.0105
		(0.0016)	(0.0016)	(0.0020)	(0.0019)	(0.0019)
	1

Table 7 :

 7 Distribution of central histology among cases and controls in the sampled strata.

			Controls			Cases
		Central histology	Subcohort	Central histology	Subcohort
	Stratum Favorable Unfavorable	fraction	Favorable Unfavorable	fraction
	1	450	2	0.27	53	4	1
	2	1,569	51	0.10	216	16	1
	4	897	17	0.13	188	20	1
				20

Table 8 :

 8 Mean results from 1,000 simultated phase-2 samples based on the NWTS data. Imputation model for unfavorable central histology: status indicator, local histology, stage, age, tumor, diameter and the interaction local histology*stage. , binary indicator of unfavorable central histology; Age0 and Age1, piecewise linear terms for age at diagnosis (years) before and after 1 year; stage, binary indicator of stage III/IV disease; diameter, diameter (cm) of excised tumor; SE, standard error.

	Weighted estimator	Standard Calibrated Re-estimated	β β (SE 4 ) (SE 4 ) β (SE 4 )	4.046 4.046 4.050	(0.537) (0.520) (0.518)	-0.669 -0.663 -0.676	(0.359) (0.324) (0.324)	0.106 0.104 0.107	(0.026) (0.017) (0.017)	1.344 1.345 1.344	(0.346) (0.273) (0.272)	0.070 0.070 0.070	(0.021) (0.015) (0.015)	-0.076 -0.076 -0.076	(0.029) (0.021) (0.021)	-2.648 -2.655 -2.651	(0.612) (0.592) (0.590)	-0.051 -0.050 -0.052	(0.051) (0.050) (0.048)
	Model 1 2 Multiple imputation Model 2 3	β (SE 4 ) β (SE 4 )	3.596 4.106	(0.555) (0.469)	-0.702 -0.537	(0.329) (0.329)	0.102 0.106	(0.017) (0.016)	1.441 1.353	(0.257) (0.251)	0.073 0.072	(0.014) (0.014)	-0.082 -0.083	(0.020) (0.020)	-2.239 -3.097	(0.611) (0.525)	-0.041 -0.065	(0.041) (0.040)
	Whole cohort 1		β	(SE)	4.042	(0.413)	-0.661	(0.326)	0.104	(0.017)	1.346	(0.244)	0.069	(0.014)	-0.076	(0.019)	-2.635	(0.464)	-0.058	(0.034)
				Parameter	UCH		Age0		Age1		Stage		Diameter		Stage*diameter		UCH*age0		UCH*age1

1

Unique estimates provided by the NWTS cohort. 2 3 Imputation model for unfavorable central histology: status indicator and local histology.

4 Asymptotic standard error (SE) of the estimate. UCH

Acknowledgments

This study was supported by a grant from the Région Île-de-France. The authors are grateful to Pierre-Yves Scarabin and Pierre Ducimetière for providing the PRIME data.


Rachid Attia

Saïd Amari
email: said.amari@lurpa.ens-cachan.fr

Claude Martinez
email: claude.martinez@irccyn.ec-nantes.fr

Control of Time-Constrained Dual-Armed Cluster Tools Using (max, +) Algebra

Keywords: Timed Event Graphs, (max, +) algebra, temporal constraints, feedback control, cluster tools

The problem studied in this paper is the control of discrete event systems subject to strict temporal constraints using (max, +) algebra. Initially we sought necessary and sufficient conditions for the existence of a causal control law guaranteeing the respect of the temporal constraints. Subsequently, a method for calculating the control law, if any, is proposed. The application which we are interested in is the control of a manufacturing semiconductor wafers process subject to strict temporal constraints.

I. INTRODUCTION:

In this work we are interested in the control, supervision, of a class of time-constrained discrete event systems modeled by deterministic Timed Event Graphs (TEG). The problem of time constraints is encountered in many industrial applications, such processes including thermal or chemical treatments [S. [START_REF] Amari | Sizing, cycle time and plant control using dioid algebra[END_REF]] , [J. [START_REF] Kim | Scheduling Analysis of Time-Constrained Dual-Armed Cluster Tools[END_REF]], the embedded systems and urban or railway transportation systems [T. [START_REF] Van Den Boom | Modeling and Control of railway networks[END_REF]]. Let us consider semi-conductor production process, once a schedule has been determined for a cluster tool, the behavior of the production system may be described as a TEG [J. [START_REF] Kim | Scheduling Analysis of Time-Constrained Dual-Armed Cluster Tools[END_REF]]. Some tasks of the production process which are executed by a robot may need to be repeated, for alignment purpose. This would perturbate the initial schedule and consequently lead to some quality loss on a part of the production. The goal of our approach is to design a control law that would be tolerant to that perturbations.

Temporal constraints have been earlier studied has a problem of verification and validation, and many authors have addressed that point of view, see for example [B. [START_REF] Berthomieu | Modeling and Verification of Time Dependent Systems Using Petri Nets[END_REF]. In our work, we consider rather, a control problem. We search for a linear feedback control law determining the firing instants of the controllable transitions, source transitions, to guarantee the respect of the temporal constraints.

In our approach the time is explicitly taken into account, and that is the main difference methods [K. [START_REF] Yamalidou | Feedback Control of Petri Nets Based on Place Invariants[END_REF]], [L. E. [START_REF] Holloway | A Survey of Petri Net Methods for Controlled Discrete Event Systems[END_REF]]. In the literature, we may also find other approaches of control using dioid algebra. In [L. [START_REF] Houssin | Control of Constrained (Max,+)-Linear Systems Minimizing Delays[END_REF] an approach based on fixed points results of antitone mappings is given, the aim of the proposed control method is to delay as less as possible the system while ensuring some given specifications. Another approach of supervision of an industrial plant is proposed in [A. M. [START_REF] Atto | Supervision of an Industrial Plant Subject to a Maximal Duration Constraint[END_REF]. In both methods, the authors consider a completely controllable TEG, i.e. all transitions are controllable, which is limiting in several real applications. Other approaches are proposed in [S. [START_REF] Amari | Control of Linear Min-plus Systems under Temporal Constraints[END_REF]] and [S. [START_REF] Amari | Design of Max-Plus Control Laws to Meet Temporal Constraints in Timed Event Graphs[END_REF]] in (min, +) and (max, +) algebras respectively. Some restrictive assumptions were considered for the control synthesis, e.g. existence of an empty path, with no tokens, from the control transition to the constrained place. The control laws were calculated under sufficient conditions.

We suggest to relax the assumptions taken by [S. [START_REF] Amari | Design of Max-Plus Control Laws to Meet Temporal Constraints in Timed Event Graphs[END_REF]] and formulate necessary and sufficient conditions for the existence of a causal linear feedback, ensuring the respect of the temporal constraints. We consider in our study the systems that can be modeled by a linear (max, +) equations subject to temporal constraints, the temporal constraints are represented by a set of (max, +) linear inequalities. We consider first the control of systems with a single control input and extend the approach to systems with multiple control inputs, assuming that the whole set of constraints is admissible. As an application, we are interested in the control of a time-constrained dualarmed cluster tool proposed in [J. [START_REF] Kim | Scheduling Analysis of Time-Constrained Dual-Armed Cluster Tools[END_REF]].

This paper is organized as follows, the second section recalls briefly the bases and tools of (max, +) algebra and TEG modeling. In the third section, we bring to light the temporal constraints problem and their formalization in (max, +) algebra. The proposed and the established results are presented. The application of the proposed method for the control of a time-constrained dual-armed cluster tool is presented in the fourth section and then we conclude.

II. PRELIMINARIES:

(max, +) algebra:

A dioid, or semiring, is a set equiped with two binary operations and called addition and multiplication, respectively. The addition is commutative, associative with identity element called ''zero''. The multiplication is also associative with identity element called ''identity', if the multiplication is commutative, the dioid is commutative'. The multiplication is distributive over addition and the ''zero'' annihilates , with respect to multiplication.

The dioid commonly called (max, +) algebra, which we consider here, is defined over the set of complete real numbers ∞, ∞ . The operations addition and multiplication correspond, respectively, to the operations maximum with ( ∞) as a "zero" element and ordinary addition with an "identity" element ( 0). The multiplication is replaced by • in the rest of this paper.

Timed Event Graphs:

Timed Event Graphs (TEGs) define a subclass of Petri nets [T. [START_REF] Murata | Petri nets, Properties, Analysis and Applications[END_REF] where each place has exactly one upstream transition and one downstream transition. i.e. this type of graphs allows modeling synchronization and parallelism phenomena but not resource sharing or mutual exclusion. An important property of TEGs is that for any circuit, the number of token is constant, therefore if all circuits are non-empty, the TEG is live. Several real systems, as flexible workshops, embedded systems and transportation systems are modeled by TEGs.

TEGs are temporized Petri nets; we denote the place relying the transition to , if any, the corresponding delay is denoted and its initial number of tokens denoted . The temporization corresponds to the minimal sojourn time of tokens in the place. The maximal marking arising in the graph is denoted . The evolution mode considered for TEGs is a maximum speed mode, i.e. a transition is fired as soon as enabled.

We define a path from a given transition to another transition as the suite of transitions and places ( , , , , ). The marking of the path is given by the sum of the number of tokens in each place of the path. The temporization of the path is also given by the sum of the path places temporizations.

Linear (max, +) model and state space representation:

It is well known that the dynamical behavior of a TEG can be expressed by a system of linear inequalities in the (max, +) algebra [F. [START_REF] Baccelli | Synchronization and Linearity[END_REF]. For this, we associate to each transition a dater, function of the integer variable , where, corresponds to the date of the firing of the transition , such that:

where: is the state variables vector, is the vector associated to the input transitions ,…, .

is the matrix whose element , equals to the temporization of the place , if any, and otherwise.

is the input matrix. Considering a maximal speed (as soon as) evolution of the TEG, the precedent inequality (I.1) can be replaced by the following equation:

The implicit equation (I.2) is usually replaced by the following explicit equation:

• • 1 I.3)
where '*' is the Kleene Star operator.

• ••• • , r times.
To get a state-space representation we must decompose each place containing more than one token to several places and transitions where each place contains at most one token. By this operation, any model (I.3) can be written in the following state-space representation:

• 1 • Matrices • and • .
III. CONTROL SYNTHESIS:

1. Temporal constraint:

The temporizations associated with places in a TEG correspond to the minimal sojourn time. In fact, the tokens are allowed to sojourn more time. For a time-constrained place a maximal sojourn time is fixed. This limitation of the maximal allowed sojourn time appears as an additional constraint that should be verified. Let us first consider a single constraint on place . We associate to this place the time interval , , is the minimal sojourn time while is the maximal one.

The constraint is expressed through the following inequality:

• , 0
where is the initial marking of the place . The inequality (II.1) is the additional constraint to be satisfied. Thus, the problem is to determine when the controllable transitions should be fired to satisfy the constraint (II.1). To consider more than a single constraint, one has to consider as many expressions (II.1) as the number of temporal constraint to satisfy.

Control synthesis:

We formulate the problem as follows; let a system given by its state-space representation (I.4) subject to the temporal constraint (II.1). Find a linear feedback • 1 , such as the constraint (II.1) being always satisfied. In order to establish necessary and sufficient conditions for the existence of a control law ensuring the respect of the temporal constraints we must discard the , trivial case where is a solution. We distinguish the two following cases:

Trivially solvable problems:

Let the autonomous system, i.e. , given by its state-space model:

• 1 hence • 1 and • 1 ,
where and are, respectively, the rows and of the matrix . Replacing and in the inequality constraint (II.1) we obtain:

• • 0 • • • 0 0 being known,
1. For a single constraint, it suffices to verify the inequality (II.2) to decide whether the problem is trivially solvable. Considering a set of q constraints, expression (II.2) has to be verified for each constraint, and furthermore, inequality (II.2.1) has to be satisfied:

• 0 • 0 , (II.2.1)
where , , H is composed of q row matrices • for l=1 to q, Q is composed of q row matrices ,

• • , and max , 1. In a recent work of [START_REF] Allamigeon | The Tropical double description method[END_REF] an algorithm is given to characterize the complete set of solutions to inequalities of the form (II.2.1). In the present work, x(0) is given.

Control synthesis for time-constrained systems:

Again, let us consider first a system with a single input and single temporal constraint, given by its (max, +) statespace model (I.4). Let being the place subject to the temporal constraint (II.1). The idea of the control is delaying the entrance of the tokens in the time-constrained place. Indeed, we block the transition by blocking the input transition . This method of control can not operate unless the blocking of the control transition does not induce an "important" blocking of the transition . The term "important" is informal, it is function of the marking and temporizations of the paths relying to and . Afterwards, we derive necessary and sufficient conditions for the existence of such control.

Let being the smallest marking of all paths relying t to , and the smallest marking of all paths from t to . The state-space representation (I.4) can be written, in the following form [S. [START_REF] Amari | Design of Max-Plus Control Laws to Meet Temporal Constraints in Timed Event Graphs[END_REF]]:

• • •
with a strictly positive integer.

Hence, for 1 the expressions of the component is given by:

• 1 • •
where is the i th row of the matrix and

• is the i th element of the vector • . As far as, for 1 is given by:

• 1 • •
Considering the definitions of (resp.

) the terms • (resp.

•

) vanish for (resp.

). Thus, we get the following expressions for and :

• 1 • • • • • 1 • •
The inequality constraint (II.1) is equivalent to

•

. Replacing and by their expressions in (II.3) we find:

The temporal constraint (II.1) is satisfied if and only if the two inequalities (II.4) and (II.5) are satisfied. From the two last inequalities (II.4) and (II.5) we derive necessary and sufficient conditions for the existence of a control law satisfying the temporal constraint (II.1) and calculate this law, if any. We must distinguish the two cases and . i) First case : In this case, we just have to verify that inequality (II.4) is satisfied. In fact, by replacing by its expression, we get:

• 1 • • • 1 • •
From which we derive the first necessary and sufficient condition, given by the following proposition:

Proposition 1: Consider a system given by its (max, +) state-space model (I.4), subject to temporal constraint (II.1). There is a linear feedback control law ensuring the constraint (II.1) if and only if the following condition holds:

• • • (II.2) (II.3) (II.4) (II.5) (II.6) • • • • 1 • • •

414

Proof: To prove the sufficiency of the proposed condition, it suffices to prove the correctness of the following implication:

• • • •
. We have:

• • • .
By multiplying the inequality by we get:

• • • • •
It suffices to choose:

• • •
• 1 , to get the equality (II.8) :

• 1 • • •
Now adding, in (max, +) sense, inequality (II.7) to equality (II.8) we get

• 1 • • • • •
. The left term, in fact, is equal to , this leads to

• • •
. Finally, we get:

• To prove the necessity of the condition (II.6) we proceed by contradiction. We suppose that the condition (II.6) doesn't hold and the constraint is satisfied and this must lead to an anomaly.

Suppose:

• • •
and:

• Inequality (II.10) leads to:

• 1 • • • •
Taking in account the expression of , two cases arise, the first is:

• • replacing in (II.11) by its expression, we get:

• 1 • • • • • • •
The second inequality of the system (II.12) contradicts the hypothesis (II.9).

The second case is:

• 1 . Replacing by its expression in (II.11), we find this time:

• 1 • • 1 • • • • 1
In this case, the first inequality of the system (II.13) leads to an absurdity. In fact, we supposed that the autonomous system doesn't satisfy the temporal constraint, see the remark in ( §III.1).

Remark1: we proved that for the case the condition (II.6) is a sufficient and necessary one for the existence of a linear feedback control law satisfying the temporal constraint (II.1). We proposed also a control law given by the following equation:

• 1 • 1 II.14)
with:

• •
ii) Second case :

In this case, in addition to inequality (II.4) treated in the first case, we deal with inequality (II.5) form which we derive the second necessary and sufficient condition. Thus, we get the same necessary and sufficient condition (II.6) and an additional condition driven from (II.5). The inequality to be satisfied is:

• • • ;
, which is the same as:

• • • ;
We have

• 1 , from which we get

• • • 1 .
Taking into account the expression of

, • 1 • •
. We derive the following condition:

• • • • • •
Proposition 2: Consider a system given by its (max, +) state-space model (I.4), subject to temporal constraint (II.1).

There is a linear feedback control law guaranteeing the respect of the temporal constraint (II.1), if and only if the two following conditions hold:

• • • (II.16) • • • • • • Proof:
The condition (II.16) arises from the case ( ), while the second condition (II.6) is a consequence of the case (

). The condition (II.17) is as a limitation of the greatest feedback . To calculate a control law we first check the condition (II.16), if it holds, we calculate a control law using the result (II.14). Then, we check condition (II.15). We don't need to solve the inequality (II.15).

The work presented here for single control input systems subject to single temporal constraint may be applied to systems that are subject to multiple temporal constraints. After derivating a control law for each constraint, assuming (II.16) and (II.17) are satisfied for each constraint, one has to check if the resulting control laws do not influence each others.

IV. APPLICATION: SCHEDULING OF TIME-CONSTRAINED DUAL-ARMED CLUSTER TOOL:

In this section we apply our control method to solve a time-constrained scheduling problem proposed in [J. Kim ]. The studied system is a semiconductor wafers manufacturing cluster tool. Processes such as some low pressure chemical vapor deposition processes require strict timing control. Unless a processed wafer leaves the chamber within a specified time limit, the wafer is subject to quality degradation due to the processing side effects. A wafer goes through different chambers and undergoes different treatments. The operations of transportation, loading and unloading to the different chambers are performed by a dual-armed handling robot. The challenge is in scheduling the operations on the handling robot to ensure the different temporal constraints of processing steps. Our approach consists in the formulation of the time-constrained scheduling problem as a control problem under strict temporal constraints. Then, we exploit the results of the proposed method to calculate a control satisfying these temporal constraints. In the following paragraphs, we give, first, a description of the manufacturing system. After, the associated TEG and (max, +) model will be given. A control law ensuring the respect of the temporal constraints is calculated, thus, a feasible schedule is found.

System description:

A cluster tool of type flow pattern (2,1) is compound of two load locks, two parallel chambers C1, C2 and a third chamber C3 in a serial configuration with C1 and C2. The same processing is performed in the chambers C1 and C2 while a different one in the chamber C3, that said, a new wafer unloaded from the load lock (LL) passes through one of the chambers C1 or C2 where a first treatment is completed then it passes to the chamber C3 where another treatment is offered. The following figure schematizes a cluster tool of type flow pattern (2, 1) and brings to light the robot work cycle.

System Modeling:

The robot work cycle is defined by the following operations sequence: unload a new wafer from the load lock move the empty arm to C1 swap the completed wafer with the unprocessed wafer at C1 move to C3 swap at C3 move to the load lock return the wafer to the load lock unload a new wafer from the load lock move to C2 swap at C2 move to C3 swap at C3 move to the load lock return the wafer to the load lock. , respectively. The linear (max, +) model of our system is given by the following linear equations:

• 1 • 1 • • • 1 • • • 1 • • • • with: 4;
1; 2; 22; 9; 1; 1.

The above system of equations can be written in a statespace representation as follows (see §II.3),

• 1

• . Matrices A and B are given below.

Control of the time-constrained system:

We mentioned above the problem of the strict time processing in the chambers. On the TEG modeling the robot work cycle, we denoted , the temporization associated to processing in chambers C1 and C2 while , for the temporization associated to the once in chamber C3. The temporal constraints are written in (max, +) algebra as follows:

ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ; 2 3
• 1 • 1 • 1 •
The analysis is performed for each input towards all of the temporal constraints by checking the conditions (II.6) for each input regarding the temporal constraints (III.1). Let being the column of the matrix , so we check the following condition:

ε ε ε ε ε • 1 ,
and are not used for control, this mean that they are vanished. To control the system we connect the TEGs resulting from and to the TEG of the system, thus we get the TEG of the controlled system.

Results analysis:

In this application, we have shown that the (max, +) method proposed in this paper can bring solutions to realworld problems encountered in industrial applications. About comparing our results to those proposed by [J. [START_REF] Kim | Scheduling Analysis of Time-Constrained Dual-Armed Cluster Tools[END_REF]], it remains not easy and less evident since the authors in [J. [START_REF] Kim | Scheduling Analysis of Time-Constrained Dual-Armed Cluster Tools[END_REF]] are interested in scheduling of dual-armed cluster tools, while we are interested in developing formal methods for control and supervision of time-constrained discrete event systems. For us, the scheduling of the dual-armed cluster tool studied here is an application between many other applications where the proposed approach can be applied.

V. CONCLUSION

In this paper we proposed a (max, +) formal method for the control of time-constrained discrete event systems. The principal contribution of this work consists in formulating necessary and sufficient conditions for the existence of linear feedback control law ensuring the respect of the strict temporal constraints. Originality lies also in the application of the proposed method for solving a time-constrained scheduling problem. We wish to generalize the method to multivariable systems, multi-inputs and several temporal constraints. Working on the application presented in this paper and many other examples, we notice that the control law found is always optimal in cycle time sense, i.e. it delays as less as possible the system while keeping the temporal constraints satisfied.

Fig 2 .

 2 Fig 2. Place subject to strict temporal constraint.

(

  The TEG of the(Fig 4) traduces the robot work cycle, each place corresponds to an operation or series of operations. We denote the temporizations associated to the time-constrained places with a closed interval of the form , where is the minimal duration of a treatment and the maximal waiting time to leave chamber after treatment ends.On the TEG the places are denoted from to and the transitions from to . Places and correspond to processing in chambers C1 and C2, respectively. While and corresponds to processing in C3. Places , , , and correspond to swap operations on which no waiting time is allowed, thus, downstream transitions of these places are uncontrollable, i.e. transitions , , and are uncontrollable. The rest of transitions are controllable, we associate to each one of them a control input. The control inputs are shown on (Fig 4). To get the (max, +) model traducing the dynamical behavior of the robot work cycle, we associate to each transition of the TEG (Fig 4) a dater denoted and to the controllable transitions , , and the control inputs , , and

Fig 3 .Fig 4 .

 34 Fig 3. Wafer flow pattern (2, 1) for dual-armed cluster tool.