Franck Corset
email: franck.corset@upmf-grenoble.fr

Laurent Doyen

Olivier Gaudoin

Bayesian Analysis of ARA Imperfect Repair Models

Keywords: Bayesian inference, Imperfect repair, Maintenance efficiency, MCMC algorithms, Repairable systems reliability, Virtual age Mathematics Subject Classification Primary 62N05, Secondary 60K10

This article proposes a Bayesian analysis of a class of imperfect repair models, the ARA models. The choice of prior distributions and the computation of posterior distributions are discussed. The presentation is unified for all ARA models and many kinds of possible priors. A numerical study on the quality of the Bayesian estimators is presented, as well as a comparison with the maximum likelihood estimators. Finally, the approach is applied to a real data set.

Introduction

The reliability of a repairable system depends on both the aging process and the efficiency of repair. The basic assumptions on repair efficiency are known as minimal repair or As Bad As Old (ABAO) and perfect repair or As Good As New (AGAN). The corresponding stochastic models for the failure process are, respectively, the Non Homogeneous Poisson Processes (NHPP) and the Renewal Processes (RP). The reality is generally between these two extreme cases: standard repair is better than minimal but not necessarily perfect. This is known as imperfect repair.

Many imperfect repair models have been proposed. The most usual are Kijima's virtual age models [START_REF] Kijima | Some results for repairable systems with general repair[END_REF], where repair is assumed to rejuvenate the system. [START_REF] Brown | Imperfect repair[END_REF] assumed that repair is perfect with probability p and minimal with probability 1p. [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF] proposed two classes of models, the ARA and ARI models, based on a reduction of virtual age or failure intensity.

Statistical inference is needed in order to estimate model parameters and compute reliability indicators from failure data. Many results exist in the NHPP and RP cases, but only a few articles deal with the statistical inference in imperfect repair models. For virtual age models, empirical studies on maximum likelihood estimators (MLE) have been presented by [START_REF] Shin | Estimating parameters of intensity function and maintenance effect for repairable unit[END_REF], [START_REF] Lim | Estimating system reliability with fully masked data under Brown-Proschan model[END_REF], [START_REF] Kaminskiy | G-renewal process as a model for statistical warranty claim prediction[END_REF], [START_REF] Gasmi | A general repair, proportional-hazards, framework to model complex repairable systems[END_REF], [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF], [START_REF] Veber | Generalized renewal process for repairable systems based on finite Weibull mixture[END_REF], [START_REF] Dijoux | A virtual age model based on a bathtub shaped initial intensity[END_REF], and [START_REF] Syamsundar | Imperfect repair proportional intensity models for maintained systems[END_REF]. Theoretical results on these estimators were derived by [START_REF] Peña | Semiparametric inference for a general class of models for recurrent events[END_REF] and [START_REF] Doyen | Asymptotic properties of imperfect repair models and estimation of repair efficiency[END_REF].

In practice, for reliable systems, only a few failures occur. So the asymptotic results on MLE cannot be used. Moreover, the engineers knowledge about the degradation and failure process could be very helpful to improve the estimations of the model parameters. Thus, a Bayesian analysis of these models is an interesting alternative to usual frequentist methods.

The Bayesian inference has been extensively studied in the NHPP case: [START_REF] Guida | Bayes inference for a non-homogeneous Poisson process with power intensity law[END_REF][START_REF] Guida | Bayes estimation of prediction intervals for a power-law process[END_REF], [START_REF] Guida | Bayesian analysis of repairable systems showing a bounded failure intensity[END_REF], [START_REF] Bar-Lev | Bayesian inference for the power-law process[END_REF], [START_REF] Campodonico | Inference and predictions from Poisson point processes incorporating expert knowledge[END_REF], [START_REF] Calabria | A reliability growth model in a Bayes decision framework[END_REF], [START_REF] Kuo | Bayesian computation for nonhomogeneous Poisson processes in software reliability[END_REF], [START_REF] Beiser | Bayes prediction for the number of failures of a repairable system[END_REF], [START_REF] Sen | Bayesian estimation and prediction of the intensity of the power-law process[END_REF], [START_REF] Ryan | Some flexible families of intensities for non-homogeneous Poisson process models and their Bayes inference[END_REF], and [START_REF] Pievatolo | Bayesian reliability analysis of complex repairable systems[END_REF]. Conversely, Bayesian analysis for imperfect repair models has been seldom studied: [START_REF] Lim | Bayesian imperfect repair model[END_REF] for the Brown-Proschan model, and [START_REF] Pan | Bayes inference for general repairable systems[END_REF] for the ARA 1 and ARI 1 models. With a different purpose, [START_REF] Sethuraman | Nonparametric Bayes estimation in repair models[END_REF] derived theoretical results for general virtual age models when the virtual ages are known.

The aim of this article is to present a Bayesian analysis of the class of ARA imperfect repair models. This class involves the ARA 1 model studied by [START_REF] Pan | Bayes inference for general repairable systems[END_REF] and also many others including the ARA model.

This article is organized as follows. The ARA models are defined in Sec. 2. The Bayesian analysis is developed in Sec. 3, with a discussion on the choice of prior distributions and the computation of posterior distributions. Section 4 studies the properties of the estimators by means of Monte-Carlo simulations. Finally, an application to real data is presented in Sec. 5.

The ARA Imperfect Repair Models

Let T i i≥1 be the successive failure times of a repairable system, starting from T 0 = 0. We assume that a repair task is performed after each failure and that repair durations are negligible. Let X i i≥1 be the times between failures X i = T i -T i-1 and N t be the number of failures observed up to time t.

Then, the failure process is a random point process. Its distribution can be characterized by the failure intensity, defined as:

∀ t ≥ 0 t = lim t→0 1 t P N t+ t --N t-= 1 t- (1)
where t-is the history of the failure process just before time t, i.e., the set of all events occurred just before t. N t-denotes the left-hand limit of N t . In most cases, the failure process is a self-excited process, i.e., t is a function of the number of failures and the failure times before t: t = t N t-T 1 T N t-. In this case, the distribution of the failure process is completely given by the intensity.

We assume that before the first failure, the intensity is a deterministic continuous function of time, denoted t , and called the initial intensity. The initial intensity characterizes the intrinsic behaviour of a new unrepaired system. is the hazard rate of the first failure time T 1 = X 1 . In the following, we assume, as in Pan Downloaded by [franck corset] at 05:32 18 September 2012 and Rigdon (2009), that the initial intensity is the intensity of a Power Law Process (PLP):

t = t -1 > 0 > 0 (2)
In general, the system is ageing intrinsically, so t is increasing and > 1.

For the two basic models, the failure intensity of a NHPP (ABAO) is a function of time:

t = t (3)
and the failure intensity of a RP (AGAN) is a function of the time elapsed since last repair:

t = t -T N t- (4)
The idea of the virtual age models proposed by [START_REF] Kijima | Some results for repairable systems with general repair[END_REF] is to assume that there exists a sequence of random variables A i i≥1 , with A 0 = 0, such that after the ith repair, the system behaves like a new one having survived without failure until A i . This property can be written:

P X i+1 > x A i X 1 X i = P Y > A i + x Y > A i A i for all x ≥ 0 (5)
where Y denotes a random variable with the same distribution as the first failure time X 1 . Then, it can easily be proved that the failure intensity is:

t = A N t-+ t -T N t- (6)
The virtual age of the system at time t is A N t-+ t -T N t-. A i is the virtual age just after the ith repair and is called the ith effective age. The NHPP is a virtual age model with ∀i A i = T i and the RP is a virtual age model with ∀i A i = 0.

This article is dedicated to the Arithmetic Reduction of Age (ARA) class of imperfect repair models, defined in [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF], which are particular virtual age models. The Arithmetic Reduction of Age model with memory m, denoted ARA m , is defined by its failure intensity:

t = t - min m-1 N t--1 j=0 1 -j T N t--j (7)
m reflects a Markovian property: it is the maximal number of previous failure times involved in the failure intensity. is a parameter which characterizes repair efficiency and will be described later. Two values of m are particularly interesting.

• m = 1. The failure intensity of the ARA 1 model is:

In this model, the effect of a repair is to reduce the virtual age just before repair of a quantity proportional to the time elapsed since last repair:

A i = A i-1 + X i -X i ⇒ A i = 1 -T i
The ARA 1 model is similar to Kijima's Type I model [START_REF] Kijima | Some results for repairable systems with general repair[END_REF] with deterministic repair effects, and is also the same as the Proportional Reduction of Age model proposed by [START_REF] Malik | Reliable preventive maintenance scheduling[END_REF]. • m = + . The failure intensity of the ARA model is:

t = t - N t--1 j=0 1 -j T N t--j (9)
In this model, the effect of a repair is to reduce the virtual age just before repair of a quantity proportional to this virtual age:

A i = 1 - A i-1 + X i
The ARA model is similar to Kijima's Type II model [START_REF] Kijima | Some results for repairable systems with general repair[END_REF] with deterministic repair effects, and is also the same as the model proposed by Brown et al. (1983).

In these models, the repair efficiency is measured by the value of . Particular values of are:

• = 1: perfect repair (AGAN); • 0 < < 1: efficient repair; • = 0: minimal repair (ABAO); and • < 0: harmful repair.

For the sake of simplicity, we will assume that repair cannot be harmful, so ∈ 0 1 . Of course, in practice, it is doubtful that repair efficiency is constant. It is possible to take into account time-dependent [START_REF] Fuqing | A general imperfect repair model considering time-dependent repair effectiveness[END_REF] or random [START_REF] Doyen | On the Brown-Proschan model when repair effects are unknown[END_REF] repair effects, but the corresponding statistical inference is much more complex. So, we focus our analysis on the rather simple ARA models, in which can be understood as a mean repair efficiency. Note that the Bayesian framework, in which a prior distribution is given on , is a way of considering non constant repair efficiency.

In this article, we derive all our results under the general framework of the ARA m model. The interest is that these results will be valid, among others, for both the ARA 1 and ARA models. For the ARA 1 model, our work extends the results obtained by [START_REF] Pan | Bayes inference for general repairable systems[END_REF]. The results for the ARA model are new.

For a point process model with parameter (here = ), the general expression of the likelihood function associated with the observation of n failures up to time t of a single repairable system is [START_REF] Cook | The Statistical Analysis of Recurrent Events[END_REF]:

L t n t n = n i=1 t i exp - n+1 i=1 t i t i-1
s ds (10) Downloaded by [franck corset] at 05:32 18 September 2012 where t n = t 1 t n denotes the vector of observed failure times, t 0 = 0 and t n+1 = t. In this situation, the process is time-truncated at time t. For a process which is failure-truncated after n failures, (10) is still valid by taking t n+1 = t n .

For the ARA m model with PLP initial intensity, the likelihood becomes, with the same kind of notations as in [START_REF] Pulcini | On the overhaul effect for repairable mechanical units: a Bayes approach[END_REF]:

L t n t n = n n V -1 exp -Z (11)
where

V = n i=1 t i -S i Z = n+1 i=1 t i -S i -t i-1 -S i S 1 = 0 and ∀i ≥ 2 S i = min m-1 i-2 j=0 1 -j t i-j-1
Let ˆ t , ˆ t , and ˆ t be the MLE of , , and , obtained by maximizing (11). It is easy to show that:

ˆ t = N t Z ˆ t ˆ t (12)
However, numerical optimization procedures are necessary to compute ˆ t and ˆ t .

Bayesian Inference

In this article, we consider a single repairable system, which is so reliable that only very few failures occur. Then, the quality of the MLE can be very poor. The Bayesian analysis improves the accuracy of parameter estimations by adding the expert knowledge to operation feedback data. In the Bayesian framework, the parameters are considered as random variables. Their prior distributions reflect the expert knowledge on the system aging and repair efficiency. The ARA models have three parameters > 0, > 0, and ∈ 0 1 . For a given prior density , the posterior density is obtained by:

n t n ∝ L t n t n (13)
The choice of prior distributions for this kind of parameters has been discussed in many articles. Non informative priors are used when no particular information is known on the failure and repair process. When such information is available, it is translated into informative priors. In practice, experts give an average behaviour and a degree of uncertainty on this information, which can be converted into prior mean and variance.

Since and are the parameters of a PLP, the prior distributions used in the Bayesian analysis of the PLP can be chosen here. Non informative priors for Downloaded by [franck corset] at 05:32 18 September 2012 and are used by [START_REF] Guida | Bayes inference for a non-homogeneous Poisson process with power intensity law[END_REF], [START_REF] Bar-Lev | Bayesian inference for the power-law process[END_REF], [START_REF] Calabria | A reliability growth model in a Bayes decision framework[END_REF][START_REF] Sen | Bayesian estimation and prediction of the intensity of the power-law process[END_REF]. Informative priors are used either by considering expert knowledge on the failure mechanism, or by using conjugate distributions. For , the gamma distribution has been used by [START_REF] Kuo | Bayesian computation for nonhomogeneous Poisson processes in software reliability[END_REF], [START_REF] Pievatolo | Bayesian reliability analysis of complex repairable systems[END_REF], and [START_REF] Pan | Bayes inference for general repairable systems[END_REF]. For , the uniform distribution has been used by [START_REF] Guida | Bayes inference for a non-homogeneous Poisson process with power intensity law[END_REF], [START_REF] Sen | Bayesian estimation and prediction of the intensity of the power-law process[END_REF], and [START_REF] Pan | Bayes inference for general repairable systems[END_REF]. The beta distribution has been used by [START_REF] Sen | Bayesian estimation and prediction of the intensity of the power-law process[END_REF].

For the repair efficiency parameter of the Brown-Proschan model, [START_REF] Lim | Bayesian imperfect repair model[END_REF] used the beta distribution. For the repair efficiency parameters of the ARA 1 and ARI 1 models, Pan and Rigdon ( 2009) also used the beta distribution.

In this article, we will choose informative, non informative, and semi informative priors on the model parameters. Informative priors on reliability indicators will also be considered.

Informative Priors on the Parameters

For the scale parameter , we choose the gamma distribution, denoted a b , with density function:

= b a a a-1 exp -b (14)
where a = + 0 t a-1 e -t dt. The expectation and variance of this distribution are, respectively, a/b and a/b2 . We will see later that this is a conjugate prior for .

For the shape parameter , we choose the uniform distribution on 1 2 , denoted 1 2 :

= 1 2 -1 1 2 (15)
is linked to the system aging in a Weibull framework. So experts are likely to give plausible values for . 1 and 2 can be understood as, respectively, the lower and upper bounds for these values. It is possible to take 1 = 1 in order to assume that, without repair, the system wears out.

For the repair efficiency parameter , since we have assumed that belongs to 0 1 , we choose the beta distribution, denoted Be c d , with density function:

= c + d c d c-1 1 -d-1 (16)
The prior expectation is c c+d and the prior variance is cd c+d 2 c+d+1 . Finally, the priors are supposed to be independent, so the joint prior density is:

= ∝ a-1 exp -b c-1 1 -d-1

Non Informative Priors

The usual non informative priors for and in the PLP case are:

∝ 1 and ∝ 1 (18)
For non harmful repair, belongs to 0 1 . So a non informative prior in this case is the uniform distribution 0 1 .

Then, with the independence of priors, the joint non informative prior is:

∝ 1 (19)

Unified Priors

Following [START_REF] Pulcini | On the overhaul effect for repairable mechanical units: a Bayes approach[END_REF], we can unify the informative and non informative cases by making the following remarks.

• For , the non informative prior is a limit case of the a b distribution, where a and b tend to zero. • For , the 1 2 density and the non informative prior can both be written:

∝ - 1 2 (20)
where = 0 corresponds to the uniform distribution, and = 1, 1 = 0, and 2 = + correspond to the non informative prior. • For , the 0 1 distribution is the same as the Be 1 1 distribution, so it is a particular case of the Be c d distribution.

So we can define a unified prior as:

∝ a-1 exp -b -c-1 1 -d-1 1 2 (21)
The informative prior ( 17) is obtained for = 0 and the non informative prior ( 19) is obtained for a = b = 0, = 1, c = d = 1, 1 = 0, and 2 = + .

It is also possible to use semi-informative priors. Indeed, the three parameters of the ARA models have well-identified meanings.

is linked to the intrinsic aging, characterizes repair efficiency, and is a scale parameter. In practice, it is much easier to obtain prior information on and than on . Therefore, one can consider, as in [START_REF] Pulcini | On the overhaul effect for repairable mechanical units: a Bayes approach[END_REF], to use the non informative prior for and the informative priors for and . This can be done by taking only a = b = 0 in (21).

Posterior Distributions

For the unified prior, the joint posterior distribution is derived using ( 13), (11), and ( 21). The posterior density of is:

n t n ∝ n+a-1 exp -b + Z n-V -1 c-1 1 -d-1 1 2
(22) Downloaded by [franck corset] at 05:32 18 September 2012

The conditional posterior densities of each parameter can also be derived as:

n t n ∝ n+a-1 exp -b + Z (23) n t n ∝ n-V -1 exp -Z 1 2 (24) n t n ∝ V -1 exp -Z c-1 1 -d-1 (25)
Equation ( 23) proves that the posterior distribution of given and is the gamma distribution n + a b + Z . Then the gamma distribution is a conjugate prior for .

Because the parameter is multidimensional and there are no conjugate priors for and , the posterior distributions will be derived by means of the Gibbs sampling algorithm, with Metropolis-Hastings steps [START_REF] Robert | Monte-Carlo Statistical Methods[END_REF]. This Markov Chain Monte Carlo (MCMC) algorithm produces samples from the joint posterior distribution, from which it is possible to compute all the features of the posterior distributions. For instance, the posterior mean, median, and mode are point Bayesian estimators of the parameters. With the posterior distribution, it is also possible to give credibility intervals for the parameters.

We propose two versions of the algorithm, one for which a prior is proposed for , and one for which is derived from (12).

Algorithm 1 (Algorithm with a Prior for ). 1. k ← 0. Choose initial values of 0 , 0 , and 0 . 2. k ← k + 1. 3. Sample k given k-1 and k-1 from (23). 4. Sample k given k and k-1 from (24). 5. Sample k given k and k from (25). 6. While k < K, go to step 2.

Since it is not possible to sample directly from ( 24) and (25), we use the Metropolis-Hastings algorithm in steps 4 and 5. For instance, in step 4, we sample from an instrumental distribution easy to simulate, with density q. Let

r = k k-1 n t n q k-1 k-1 k k-1 n t n q k-1 (26)
Then, k will be equal to with probability min r 1 and to k-1 otherwise. For a symmetric transition, q k-1 = q k-1 , so (26) becomes

r = k k-1 n t n k-1 k k-1 n t n (27)
Here, we take for q a Gaussian distribution centered on the given value of .

As said before, it is not easy to give a prior on the scale parameter . A first possibility is to choose a non informative prior by taking a = b = 0 in (21). Another possibility is to use the fact that, in (12), the MLE of is computed as a function of the MLE of and . Then, we propose a second version of the algorithm, for which priors are given only for and , and is computed as a function of and , according to (12).

Algorithm 2 (Algorithm with no Prior for ). 1. k ← 0. Choose initial values of 0 and 0 . Downloaded by [franck corset] at 05:32 18 September 2012 2. Compute 0 as a function of 0 and 0 from (12). 3. k ← k + 1. 4. Sample k given k-1 and k-1 from (24). 5. Sample k given k-1 and k from (25). 6. Compute k as a function of k and k from (12). 7. While k < K, go to step 3.

Informative Prior on Reliability Indicators

Giving a prior is facilitated if the considered quantity has a clear physical meaning, such as a reliability indicator. In the PLP case, [START_REF] Guida | Bayes inference for a non-homogeneous Poisson process with power intensity law[END_REF], [START_REF] Calabria | A reliability growth model in a Bayes decision framework[END_REF], and [START_REF] Sen | Bayesian estimation and prediction of the intensity of the power-law process[END_REF] proposed a gamma prior for the expected number of failures E N t . With an informative prior on , it is possible to derive a prior for . For the ARA model, the mean number of failures takes into account both the intrinsic aging (through and ) and repair efficiency (through ). Since these parameters can have compensating effects, it is not so easy to give a prior value for E N t . Then, we propose to give a prior on an indicator which is linked only to the intrinsic aging, the mean time to the first failure:

= E T 1 = 1 1/ 1 + 1 ( 28
)
Since appears at the denominator of , we choose an inverse gamma prior, denoted I A B , for . Then, the prior density is:

= B A A -A-1 exp - B ( 29
)
For A > 2, the expectation and variance of this distribution are, respectively, B

A-1 and

B 2 A-1 2 A-2 .
The conditional prior density of is derived thanks to a change of variables:

= B A A/ -1 A 1 + 1 A exp - B 1/ 1 + 1 (30)
If we keep the informative priors 1 2 for and Be c d for , the joint prior distribution is:

= ∝ A/ -1 1 + 1 A exp - B 1/ 1 + 1 c-1 1 -d-1 1 2 (31)
and the joint posterior distribution is:

n t n ∝ n+A/ -1 1 + 1 A n-1 c-1 1 -d-1 V -1 exp -Z - B 1/ 1 + 1 1 2
(32) Downloaded by [franck corset] at 05:32 18 September 2012

Posterior Distributions of Reliability Indicators

With the posterior distribution of the parameters, it is possible to determine the posterior distribution of any quantity of interest, function of these parameters, such as reliability indicators. First, we give the posterior distribution of the mean time to the first failure = E T 1 studied previously. Then, following [START_REF] Pulcini | On the overhaul effect for repairable mechanical units: a Bayes approach[END_REF], we derive the posterior distribution of the failure intensity and the expected number of failures. In this section, we do not consider the case where is computed as a function of and .

3.6.1. Posterior distribution of the mean time to the first failure. For the estimation of = E T 1 we consider two cases. In the first one, we use the prior distributions on the parameters, and in the second one we use the prior distribution on .

In the first case, from Eqs. ( 22) and ( 28) with a change of variables, the posterior distribution of is found to be:

n t n ∝ 2 1 1 0 n-+1 V -1 1 + 1/ n+a c-1 1 -d-1 -n+a -1 exp -b + Z 1 + 1/ d d (33)
Thus, the Bayesian estimator of is the posterior mean:

E n t n = 1 C 2 1 1 0 n-V -1 1 + 1/ c-1 1 -d-1 n + a -1/ b + Z n+a-1/ d d (34)
where

C = 2 1 1 0 n-V -1 c-1 1 -d-1 n + a b + Z n+a d d
The integrals involved in these equations are computed by Monte Carlo simulations. It is also possible to obtain credibility intervals for from Eq. ( 33).

Secondly, we consider the case of Sec. 3.5, where the inverse gamma prior ( 29) is chosen for . With the informative priors 1 2 for and Be c d for , the posterior distribution of is:

n t n ∝ 2 1 1 0 n V -1 1 + 1/ n c-1 1 -d-1 -n-A-1 exp -Z 1 + 1/ - B d d (35)
The Bayesian estimator of is the posterior mean. Downloaded by [franck corset] at 05:32 18 September 2012 3.6.2. Posterior distribution of the failure intensity. The failure intensity at a time

∈ t i t i+1 is = -S i+1 -1
A change of variables in ( 22) leads to the posterior distribution of :

n t n ∝ 2 1 1 0 V -1 c-1 1 -d-1 a+ -S i+1 n+a -1 n+a-1 exp - b + Z -S i+1 -1 d d (36)
The Bayesian point estimator of is the posterior mean:

E n t n = 1 C 2 1 1 0 V -1 c-1 1 -d-1 n-+1 -S i+1 -1 b + Z n+a+1 d d ( 37
)
where C is a normalizing constant. The integrals involved in (37) are computed by Monte Carlo simulations. Credibility intervals for can also be obtained from (36). Then, computing these values for each , it is possible to draw a plot of the Bayesian point and interval estimators of the failure intensity.

Posterior distribution of the expected number of failures.

A usual result of the theory of point processes is that the expected number of failures at any time is linked to the failure intensity by:

N = 0 s ds + M (38)
where M ≥0 is a martingale. Having observed failure times t 1 < • • • < t n , for ∈ t i t i+1 , 1 ≤ i ≤ n, the cumulative failure intensity is:

= 0 s ds = i j=1 t j t j-1 s ds + t i s ds (39)
For the ARA models, we obtain:

= Z i (40)
where

Z i = i j=1 t j -S j -t j-1 -S j + -S i+1 -t i -S i+1 (41)
Of course, Z n t = Z . Then, it is possible to make a Bayesian prediction of the expected number of failures at time by computing the posterior distribution of . For a given data Downloaded by [franck corset] at 05:32 18 September 2012 set, comparing the prediction with the observed value of N will give an indication of the goodness-of-fit of the ARA model.

Note that the expression of is similar to that of = -S i+1 -1 . In fact, the results for can be derived from those on by replacing -S i+1 -1 by Z i . Then, the posterior distribution of is:

n t n ∝ 2 1 1 0 V -1 c-1 1 -d-1 n- Z i n+a n+a-1 exp - b + Z Z i d d (42)
and the Bayesian point estimator of is the posterior mean:

E n t n = 1 C 2 1 1 0 V -1 c-1 1 -d-1 n-Z i b + Z n+a+1 d d ( 43
)
where C is a normalizing constant. For = t and i = n, Z n t = Z . With the non informative prior for (a = b = 0), it is easy to check that E t n t n = n.

Simulation Results

Comparison of MLE and Bayesian Estimators

Here, we compare the maximum likelihood and Bayesian estimators by means of Monte Carlo simulations. We simulate 200 samples of size n ∈ 5 10 25 of both ARA 1 and ARA models, with parameters equal to = 0 01, = 3, and = 0 6. These values are chosen in order to reflect the behaviour of a system with important aging ( = 3), good repair efficiency ( = 0 6) and a mean time of the first failure around 4 (E T 1 = 4 14).

For each sample, we compute the MLE and several Bayesian estimators of each parameter, corresponding to different choices of priors. We also compute the estimators of E T 1 .

The prior for is non informative ( = 1 1 = 0 2 = + ) or uniform on 1 4 ( = 0 1 = 1 2 = 4). This is a rather good prior since the prior mean is 2.5 and the prior standard deviation is 0.87.

Three priors are chosen for : the non informative (c = d = 1), a good one (c = 1 652 d = 0 708 with a prior mean equal to 0.7), and a bad one (c = 0 708 d = 1 652 with a prior mean equal to 0.3). For both informative priors, the standard deviation is 0.25.

For , we consider several possibilities:

• The non informative prior: a = b = 0.

• A good informative prior: a = 2 25 b = 150, with a prior mean equal to 0.015 and a prior standard deviation equal to 0.01. • The estimator of computed as a function of the estimators of and , according to (12). • The prior for given computed from a prior on E T 1 , according to (30).

The prior for E T 1 is inverse gamma with A = 8 25 and B = 36 25, leading to a prior mean equal to 5 and a prior standard deviation equal to 2. Downloaded by [franck corset] at 05:32 18 September 2012 In these figures, the following notations are adopted. The corresponding values of the hyperparameters are given in Table 1.

• MLE: Maximum Likelihood Estimator.

• NI: Non Informative prior on each parameter.

• GP: Good Prior for , prior for .

• GC: Good prior for , Computed.

• GNI: Good prior for , non Informative for .

• GT1: Good prior for , prior on E T 1 .

• BP: Bad Prior for , prior for .

• BC: Bad prior for , Computed.

• BNI: Bad prior for , Non Informative for .

• BT1: Bad prior for , prior on E T 1 .

The behavior of the estimators is essentially similar for the ARA 1 and ARA models. The decrease of the MSE is not always true for so small samples. The worst estimators are the non informative and the MLE. For both estimators of , the mse is so large that it does not appear in Figs. 3 and4. Both estimators overestimate strongly , but the MLE is correct for and , while the non informative is correct Downloaded by [franck corset] at 05:32 18 September 2012 only for . As expected, the estimators are globally better when the prior on is good (G) than when it is bad (B).

The comparison of the ways to estimate leads to the logic conclusion that the best estimators are obtained for a good informative prior (P). The three other estimators have similar behaviors, which slightly differ depending on the model and parameters. The estimator with a prior on E T 1 (T1) can be used if a good knowledge on this value is available. The estimator with a non informative prior on (NI) and the estimator with computed (C) do not need any kind of knowledge on . Both need the MCMC algorithm. C is simpler to compute than NI because the MCMC algorithm involves only two parameters instead of three. Downloaded by [franck corset] at 05:32 18 September 2012 For the estimation of in the ARA 1 model, the prior has a strong impact. We can see in Fig. 5 that the mse of all the estimators with good prior are much lower than the mse of all the estimators with bad prior. For the ARA model, this result is still true for small n, but much weaker.

For the estimation of E T 1 , it is interesting to note that the best estimator is logically obtained when a prior on E T 1 is used (GT1). The worst estimator is the MLE, followed by NI. Here again, the G estimators are better than the B estimators.

Table 2 gives the relative mean squared error for the model parameters and for E T 1 . For a quantity , the relative MSE is MSE/ 2 , where mse is taken as an average value of the MSE previously computed. The less this value is, the better is estimated. From this table, we see that is much less well estimated than the other parameters. But it doesn't prevent E T 1 from being well estimated.

Finally, the Bayesian approach will be particularly profitable for the estimation of and E T 1 . Of course, the greater is the knowledge on the parameters, the better are the estimations. But Bayesian estimators are still satisfactory even with a bad prior.

Study of One Data Set

In practice, one has to analyze one data set made of n successive failure times of a repairable system. The first objective is to estimate the parameters of ARA models and reliability indicators. Thanks to the Bayesian analysis, it is also possible to compute credibility intervals on these quantities (the equivalent is not possible with the MLE). We can also compare prior and posterior distributions. In this section, we apply this approach to data simulated from the ARA 1 model. In the next section, it will be applied to a real data set for both ARA 1 and ARA models.

Table 3 presents n = 5 failure times simulated from the ARA 1 model with parameters = 0 01, = 3, and = 0 6.

Tables 4-7 present point estimates of the parameters , , , and of E T 1 . The considered estimators are the MLE and the Bayesian GP, GNI, BNI. For the Bayesian estimators, the 80% credibility intervals and posterior standard deviations are also given.

The MLE is naturally of bad quality since we have to estimate three parameters with five observations. Indeed, the MLE of , , and E T 1 are very far from the These figures show that the Bayesian estimate provides a better fit than the MLE: the blue line is closer to the red line than the green line. The real intensity, cumulative intensity, and observed number of failures are always within the credibility limits. In Fig. 12, the ML and GNI estimates are equal to n for t = t n . For GNI, this property has been stated in Sec. 3. For the MLE, this property comes from ( 12) and (40).

Application to Automobile Data

In this section, we present an application to automobile data studied by [START_REF] Ahn | Estimating parameters of the power-law process with two measures of time[END_REF] and [START_REF] Guida | Bayesian analysis of repairable systems showing a bounded failure intensity[END_REF]. The data given in Table 8 are n = 18 failure times of an AMC Ambassador car owned by the Ohio state government. The MLEs of the parameters for the PLP, ARA 1 , and ARA models are given in Table 9. Usual model selection criteria are also given. For a model with k estimated parameters and n observed data, these criteria are defined as:

• max ln is the maximal log-likelihood;

• AIC = -2 max ln + 2k;

• AICc = -2 max ln + 2kn n-k-1 ; and • BIC = -2 max ln + 2k ln n.

The "best" model has the smallest AIC/AICc/BIC. For all the criteria, the ARA 1 model is found to be better than the PLP. For all the criteria except the BIC, the ARA model is also better than the PLP. The best model is ARA 1 . The estimates of are close to 3, indicating a strong intrinsic wear-out. The estimates of and in the two models are very different. It is not surprising since the models have different meanings (considerations of that kind are given in [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF]. Anyway, the estimates of indicate a good repair efficiency. Now we perform the Bayesian analysis of these data. Since we have no prior knowledge on this system, the first idea is to use the completely non informative Bayesian estimators. But it makes sense to assume that this automobile is wearing. So we will also use the uniform prior on 1 4 for . For , it is interesting to note that the MLE are relatively close to some values used as prior means in the simulation study: for the ARA 1 model, ˆ t = 0 898 is not far from 0.7 (the GNI of Sec. 4.1) and for the ARA model, ˆ t = 0 246 is close to 0.3 (the BNI of Sec. 4.1). Then, we will also use the corresponding informative priors for .

Tables 11-16 present point estimates, the 80% credibility intervals and posterior standard deviations of the parameters , , and for ARA 1 and ARA models. In these tables, the following estimates are used. The corresponding values of the hyperparameters are given in Table 10.

• MLE: Maximum Likelihood Estimator.

• NI: Non Informative prior on each parameter.

• NIU: Non Informative prior for and , uniform prior for .

• GNI: Non Informative prior for , uniform prior for , beta prior for for the ARA 1 model. • BNI: Non Informative prior for , uniform prior for , beta prior for for the ARA model.

For the ARA 1 model, the NI estimates are very different from the others, especially for . They also have rather high standard deviation. The GNI estimates are globally closer to the MLE than the NIU, with smaller standard deviation and narrower credibility intervals. This is due to the fact that we have chosen a prior close to the MLE. For the ARA model, these properties are no longer true. The NIU estimates are closer to the MLE than the BNI. Downloaded by [franck corset] at 05:32 18 September 2012 The NIU estimates give coherent results and are not submitted to the subjective choice of a prior for , so we recommend them in this case.

Figures 13 and 14 present, respectively, the NIU Bayesian estimates of the intensity and expected number of failures, for the ARA 1 model. As before, the Bayesian point estimate is in solid line, the credibility limits are in (blue)/dashed lines and the MLE is in (blue)/dotted (green) line.

In these figures, the Bayesian estimates are very close to the ML estimates, even if the estimates of the parameters are not close. In fact, this result is also true with the NI and GNI estimates. Moreover, we can see in Fig. 14 that the estimates are also very close to the observed number of failures. On the one hand, this means that all the estimation methods manage to catch the behavior of the failure process and to estimate correctly reliability indicators. On the other hand, the ARA 1 model provides a good fit to these data.

Conclusion and Future Work

In this article, we presented a unified approach of the Bayesian analysis of ARA imperfect repair models, including all ARA models and many kinds of possible priors. It allows to compute meaningful reliability indicators (with credibility intervals) even when very few data are available. Other ideas can be explored in the same framework: computing other reliability indicators such as predictive MTTF, performing Bayesian tests, analyzing several data sets in parallel, etc.

Another promising area is the Bayesian inference of imperfect maintenance models, where both preventive and corrective maintenance are performed. Until now, this case has been studied by [START_REF] Pulcini | On the overhaul effect for repairable mechanical units: a Bayes approach[END_REF] and [START_REF] Sheu | A Bayesian approach to an adaptative preventive maintenance model[END_REF], with only ABAO CM effects, and by [START_REF] Yu | Parameter estimation for a repairable system under imperfect maintenance[END_REF] with ARA 1 CM and PM effects.

Figure 1 .

 1 Figure 1. ARA 1 -mse for .

Figure 2 .

 2 Figure 2. ARA -mse for .

Figure 3 .

 3 Figure 3. ARA 1 -mse for .

Figure 4 .

 4 Figure 4. ARA -mse for .

Figure 5 .Figure 6 .

 56 Figure 5. ARA 1 -mse for .

Figure 7 .Figure 8 .

 78 Figure 7. ARA 1 -mse for E T 1 .

  compare the prior density (solid line) and the posterior density (histogram from the MCMC samples) of the GNI and BNI estimators of . In both cases, the posterior is brought closer to the MLE. The posterior distribution is

Figure 9 .

 9 Figure 9. GNI-Comparison of prior and posterior densities. (color figure available online.)naturally more concentrated around the true value = 0 6 for the good prior than for the bad prior.Finally, the GNI Bayesian estimates of the failure intensity and expected number of failures are plotted in Figs. 11 and 12. The solid line represents the Bayesian point estimate, the credibility limits are in dashed (blue) lines. The true value and maximum likelihood estimates of the functions are, respectively, plotted in dash-dot (red) line and dotted (green).These figures show that the Bayesian estimate provides a better fit than the MLE: the blue line is closer to the red line than the green line. The real intensity, cumulative intensity, and observed number of failures are always within the credibility limits. In Fig.12, the ML and GNI estimates are equal to n for t = t n . For GNI, this property has been stated in Sec. 3. For the MLE, this property comes from (12) and (40).

Figure 10 .

 10 Figure 10. BNI-Comparison of prior and posterior densities. (color figure available online.)

Figure 13 .

 13 Figure 13. AMC data-ARA 1 -NIU estimate of the intensity. (color figure available online.)

Figure 14 .

 14 Figure 14. AMC data-ARA 1 -NIU estimate of the expected number of failures. (color figure available online.)

Table 1

 1 Hyperparameter values for the Bayesian estimates

		a	b		1	2	c	d	A	B
	NI	0	0	1	0	+	1	1
	GP	2.25	150	0	1	4	1.652	0.708
	GC			0	1	4	1.652	0.708
	GNI	0	0	0	1	4	1.652	0.708
	GT1			0	1	4	1.652	0.708	8.25	36.25
	BP	2.25	150	0	1	4	0.708	1.652
	BC			0	1	4	0.708	1.652
	BNI	0	0	0	1	4	0.708	1.652
	BT1			0	1	4	0.708	1.652	8.25	36.25

Table 2

 2 Relative mean squared error

	E T 1

Table 3

 3 Simulated data set from the ARA 1 model

	Downloaded by [franck corset] at 05:32 18 September 2012
	6.17	8.30	9.20	10.47	12.84

Table 4

 4 One data set-ARA 1 -estimation of but the MLE of is quite good. The GP estimates (good priors for and ) are always good, in mean and variance, as well as for the length of credibility intervals. The GNI is a good estimate for and . The BNI is a good estimate for and . For , the BNI is highly biased, because it is close to the (bad) prior mean. Figures 9 and 10

	Point estimate	Lower limit	Upper limit	Standard deviation
	MLE	0.00003
	GP	0.0123	0.00360	0.0236	0.0086
	GNI	0.0277	0.00095	0.0718	0.0648
	BNI	0.0122	0.00053	0.0357	0.0270
			Table 5
		One data set-ARA 1 -estimation of
	Point estimate	Lower limit	Upper limit	Standard deviation
	MLE	5.869
	GP	2.921	2.344	3.572	0.465
	GNI	2.954	1.826	3.825	0.755
	BNI	3.070	2.030	3.829	0.655
	true values,

Table 6

 6 One data set-ARA 1 -estimation of

	Downloaded by [franck corset] at 05:32 18 September 2012	Point estimate	Lower limit	Upper limit	Standard deviation
	MLE	0.554
	GP	0.605	0.339	0.821	0.192
	GNI	0.538	0.242	0.818	0.214
	BNI	0.317	0.040	0.629	0.215
			Table 7
		One data set-ARA 1 -estimation of E T 1
		Point estimate	Lower limit	Upper limit	Standard deviation
	MLE	5.424
	GP	4.619	3.5939	5.8197	0.9524
	GNI	4.7645	3.2231	6.4155	1.2888
	BNI	5.6575	3.8603	7.5476	1.4780

Table 9

 9 AMC data-Maximum likelihood estimates and model selection criteria

		ˆ t	ˆ	t	ˆ t	max ln	AIC	AICc	BIC
	PLP	1 32 10 -4	1.625		-95.15	194.29	195.09	201.86
	ARA 1 ARA	1 30 10 -7 2 12 10 -9	3.102 3.583	0.898 0.246	-92.00 -92.68	189.99 191.36	191.71 193.07	201.33 202.70

Table 10

 10 Hyperparameter values for the Bayesian estimates

		a	b		1	2	c	d
	NI	0	0	1	0	+	1	1
	NIU	0	0	0	1	4	1	1
	GNI	0	0	0	1	4	1.652	0.708
	BNI	0	0	0	1	4	0.708	1.652
				Table 11
			AMC data-ARA 1 -estimation of
		Point estimate	Lower limit		Upper limit	Standard deviation
	MLE	1 30 10 -7
	NI	4 87 10 -4		1 20 10 -5		1 44 10 -3	7 07 10 -4
	NIU	1 03 10 -5		6 09 10 -7		3 19 10 -5	1 99 10 -5
	GNI	3 05 10 -7		1 21 10 -9		6 73 10 -7	2 78 10 -7

  Downloaded by [franck corset] at 05:32 18 September 2012

			Table 12
		AMC data-ARA 1 -estimation of
		Point estimate	Lower limit	Upper limit	Standard deviation
	MLE	3.102
	NI	1.736	1.362	2.184	0.331
	NIU	2.430	1.989	2.696	0.254
	GNI	3.104	2.718	3.894	0.418
			Table 13
		AMC data-ARA 1 -estimation of
		Point estimate	Lower limit	Upper limit	Standard deviation
	MLE	0.898
	NI	0.680	0.284	0.939	0.256
	NIU	0.801	0.626	0.938	0.159
	GNI	0.866	0.800	0.928	0.054
			Table 14
		AMC data-ARA -estimation of
		Point estimate	Lower limit	Upper limit	Standard deviation
	MLE	2 12 10 -9
	NI	2 60 10 -6	2 00 10 -7	8 20 10 -6	4 41 10 -6
	NIU	2 73 10 -7	5 41 10 -9	8 41 10 -7	6 79 10 -7
	BNI	4 12 10 -5	4 83 10 -6	1 09 10 -4	4 79 10 -5

t = t -T N t- (8) Downloaded by [franck corset] at 05:32 18 September 2012

(17) Downloaded by [franck corset] at 05:32 18 September 2012

Downloaded by [franck corset] at 05:32 18 September 2012

Acknowledgment

This research is supported by a grant from EDF R&D. The authors thank Emmanuel Remy for his helpful comments.


Charlotte Nivon

Luc Fontana

Béatrice Trombert-Paviot

Jean Noel Bally

Dr Carole Pelissier
email: pelissier.carole.chu@gmail.com.ing

General practitioners' management of patients consulting for "mental distress at work"

Keywords: Medical consultations, psychosocial risk factors, work-related common mental disorders

Introduction

Among work-related diseases, mental health problems are the second most frequent group after musculoskeletal disorders and the first cause of workrelated sickness [START_REF] Valenty | Surveillance Programme of Work-related Diseases[END_REF][START_REF] Ikonen | Work-Related Primary Care in Occupa-444 tional Health Physician's Practice[END_REF]. Previous studies have shown that common mental disorders (CMD) can be related to work characteristics such as: job insecurity, work-Rivière et al., in primary care, prevalence estimates of CMD range from 3% to 25% for anxiety disorders and 6% to 25% for depression [START_REF] Rivière | Which work-related characteristics are most strongly associated with common mental disorders? A cross-sectional study[END_REF]. General practitioners (GPs) are key players in the management of patients consulting for mental distress at work, setting up individualized medical follow-up with specialist referral if need be. On the frontline for individuals with mental health problems, they often deal with work-related common psychiatric disorders [START_REF] Rivière | Prevalence of work-related common psychiatric disorders in primary care: The French Héraclès study[END_REF]. To our knowledge, a few studies have evaluated the prevalence of workrelated mental disorders [START_REF] Rivière | Prevalence of work-related common psychiatric disorders in primary care: The French Héraclès study[END_REF][START_REF] Hussey | Work-related ill health in general practice, as reported to a UK-wide surveillance scheme[END_REF].

The present study assessed the frequency of consultation for mental distress at work, GPs' medical management and the factors underlying consultation.

Methods

The study was approved by the French data protection commission on November 21, 2016, and by University Hospital ethics committee on January 13,

(IRBN032017/CHUSTE). A self-administered

questionnaire was e-mailed to all GPs (n = 800) in the Loire administrative department of France via the local Order of Medicine (OM), to be filled out on-line using the LimeSurvey application. Data were secured and anonymized. Respondents could make free comments at the end of the questionnaire. GPs were contacted for this study, and were asked to include patients on the following criteria:

-male or female aged > 18 years; -first consultation for mental distress at work.

Non-working patients and those unable to read

French were excluded.

GPs received a document explaining the study procedure, inclusion period and inclusion and exclusion criteria, and were asked to provide patients with clear, honest and appropriate oral information; eli-95 gible patients also received a written information 96 document. Patients were free to decline to take part 97 in the study without impact on treatment.

98

The patient questionnaire covered: or a psychiatrist (Table 3). In regard to mental distress at work, burnout is not a disease recognized in the international classifications (ICD-10 and DSM-5) or in the official French list of occupational diseases. Regional Occupational Disease Recognition Committees (Comités Régionaux de Reconnaissance des Maladies Professionnelles:

CRRMP) may make decisions on a case-by-case basis, but patients have to show at least 25% permanent incapacity in order to be examined. There

  Survey variables comprised sociodemographic data (age, gender) and occupational data (working time per week, number of consultations for mental distress at work, type of medical treatment). The rate of consultations for mental distress at work was determined by dividing the number of patients consulting the GP per week by the number examined for mental distress at work. A transversal descriptive study was made of patients in the Loire area consulting their GP for mental distress at work from February to end May 2017, using two self-administered paper questionnaires: one survey was completed by the patient and one survey was completed by the GP. Twenty-two

167 4 .

 4 Discussion168The frequency of consultations for mental distress 169 at work was about 2%. Although a majority reported 170 impact on family life and health, patients may wait 171 several weeks or months before consulting.172The rate of GP consultation for mental distress at 173 work was about 2%. Héraclès[START_REF] Rivière | Which work-related characteristics are most strongly associated with common mental disorders? A cross-sectional study[END_REF], in a series of GPs 174 in the Nord-Pas-de-Calais Region of France, between of work-related psychological disorder in general 177 practice: a quarter of working patients presented 178 work-related psychological suffering. Analysis of 179 treatment testified to strong GP commitment, with 180 psychological interviews and medical prescriptions 181 in most cases and sick-leave in half, but with little 182 referral to occupational physicians, psychologist or 183 psychiatrists. The present results likewise found little 184 communication between GPs/patients and occupa-185 tional physicians regarding mental distress at work. 186 A health and job retention policy based on a worker 187 support approach should include practices targeting 188 primary and secondary prevention [9]. According to 189 the French Labor Code (article L. 4622-3), occupa-190 tional physicians play a purely preventive role, to 191 avoid any impairment of workers' health status due 192 to working conditions, notably by monitoring health 193 and safety conditions at work, contagion risk and 194 health status. The occupational physician acts as a 195 consultant for the management, employees, staff rep-196 resentatives and the health and safety and working 197 conditions committee (Comité d'Hygiène de Sécurité et des Conditions de Travail), providing advice on improving living and working conditions in the firm, adapting work-posts, techniques and work rhythm to actual human physiology, protecting workers against nuisances in general and especially work accident risks or use of dangerous substances. The French health authority (Haute Autorité de Santé), in a 2017 report entitled "Clinical identification and management of burnout syndrome" [10], recommends that, with the patient's agreement, the GP liaise with the occupational physician or an occupational pathology specialist, to raise an alert and have some information about the workplace. Analysis of the work-post and working conditions by a multidisciplinary team coordinated by the occupational physician is indispensable. Factors hindering cooperation between GPs and occupational physicians include medical confidentiality and concerns about the patient's job security [11].The individual medical records are not presently available to the occupational physician, and exchanges with the GP require the employee's agreement. Few GPs refer their patients to occupational physicians during sick-leave for the purposes of the end-of-leave consultation [12]. To promote return to work after sick-leave exceeding 3 months, a pre-return consultation with the occupational physician is set up by the GP, the national insurance consultant physician or the patient (Labor Code, Art. R.4624-29), in which the occupational physician may recommend workpost adaptation or change or occupational training to facilitate change of work-post or job. Unless the employee objects, the occupational physician informs the employer and the consultant physician of these recommendations so that all necessary measures should be taken to ensure return to work (Labor Code, Art. R. 4624-31). According to the present findings, few patients spoke about their mental distress at work to any member of the occupational health department before consulting their GP. There are numerous factors hindering communication between employees and occupational physicians: misapprehension of the role of the occupational physician and of the possibility of meeting, difficulty of making contact, lack of a relation of confidence, and concerns about independence with respect to the employer [12]. Employees often mistrust occupational physicians due to the influence they have over hiring decisions through aptitude assessment and their role as advisors to the employer [12]. Occupational physicians, however, actually play 250 a key role in keeping the employee in work. 251 According to the Labor Code (Art. R. 4624-31 252 and R. 4624-32), return-to-work consultations with 253 an occupational physician are set up after maternity 254 leave, sick-leave for occupational disease, or more 255 than 30 day's leave for a work accident or non-256 occupational illness or accident. 257 The return-to-work examination aims: 258 -to check that the original or revised work-post is 259 compatible with health status; 260 -to examine work-post alterations or adapta-261 tions or job change proposed by the employer 262 in the light of any occupational medicine rec-263 ommendations made following the pre-return 264 consultation; 265 -to recommend work-post alterations or adapta-266 tions or job change; 267 -or, as appropriate, to recommend the employee 268 as inapt to return to work. 269 The issue of mental distress at work highlights 270 the role of psychosocial risks in the deterioration 271 of employees' mental health [13]. More than three-272 quarters of GP respondents to a web questionnaire 273 thought GP training ion mental distress at work to 274 be inadequate: "Medical school culture does not deal 275 with issues of health at work" [14]. Notably, no work 276 accidents were reported in the present sample, despite 277 patients often reporting trigger factors leading to 278 GP consultation. Certain recent misadventures with 279 the Order of Medicine, which transmitted employ-280 ers' complaints against practitioners who made the 281 connection between work and ill health, may cause 282 physicians to hesitate to make such work-accident 283 reports. In this context, support from the occupational 284 physician or national health insurance consultant is 285 precious [14]. In work accidents, the relation with 286 working conditions is implicit, and the health insur-287 ance authorities do not expect the GP to demonstrate 288 it. The Regional Union of Community Health Pro-289 fessionals (Union Régionale des Professionnels de 290 Santé Médecins Libéraux) of the Provence-Alpes-291 Côte-d'Azur Region of France drew up a guide to 292 distress at work: "Le médecin libéral face à la souf-293 france au travail de ses patients" ("The community 294 physician and patients' distress at work"), listing mis-295 takes to be avoided concerning the terms employed in 296 sick-leave or work-accident certificates: the medical 297 certificate must be strictly descriptive and clinical, 298 dealing with the pathology rather than its cause; 299 thus, the physician should not write anything like 300 "depressive reaction to working conditions". Equally, "harassment" is a legal term, not a medical diagnosis. No declarations of occupational disease were made in the present series. In a 2007 study, "Les médecins face à la santé au travail : une étude qualitative dans le Sud-Est de la France" ("Physicians and health at work: a qualitative study in South-East France"), under-declaration of occupational disease concerned occupational cancer, asthma and musculoskeletal disorder [15]. Respondent physicians stressed the difficulty of drawing up the initial medical certificate.

  are no occupational disease tables in France relating psychosocial risk and mental pathologies; however, a major depressive episode or anxiety disorder may be recognized as an occupational disease by the CRRMP in an "exceptional" procedure provided for in §4, Article L 461-1 of the Health Insurance Code. Patients' expectations consist above all in wanting their doctor to hear them out. Listening is a fundamental pillar of treatment. According to a 2017 article in the Revue Médicale de Liège, "Professional fatigue syndrome (Burnout)", only attentive listening allows the different clinical presentations, somatic and psychological, of distress at work to be interpreted [16]; the authors recommend scheduling a second consultation, to leave the patient more time. More than 80% of physicians in the present study reported that consultations involving distress at work are longer than usual. A large-scale quantitative and qualitative study conducted by the Technologia consultancy under the auspices of the UMC health insurance company in 2012, entitled "Les effets du travail sur la vie privée" (Impact of work on private life) sought to understand and, if possible, measure the impact of working life on private life, in all dimensions: social life, family life, love life, and sex life. Work monopolized daily life for 50% of employees (and 61% of executives), ahead of family life (30%), the life of the couple (17%) and social life (3%). According to Technologia, most employees are torn between the importance they give to their work and the importance they struggle to give to their private life, and family life in particular: working time tends to creep up on the rest. In the present study, most patients reported history of distress going back several months or even years. By initiating more systematic dialogue about their patients' work, GPs could help improve prevention of occupational risk, in teamwork with the occupational physician. Most patients in the study were expecting sickleave, which all were accorded, for less than two weeks in most cases. It seems to be essential to get away from work and its anxiety-inducing environment so as to be able to take some distance. Sick-leave duration varies with the severity of distress. A burnout prevention guide was drawn up by researchers and field experts in collaboration with the INRS (Institut National de Recherche et de Sécurité: national research and safety institute) and ANACT (Agence Nationale pour l'Amélioration des Conditions de Travail: national agency for improvement in working conditions), advising 2-3 months' sick-leave associated to psychotherapy. In reality, sick-leave averages between 6 and 18 months, according to Marie Peze, a Doctor of Psychology, psychoanalyst and expert witness who set up the first "distress at work" consultation in the Nanterre hospital reception and care center in 1997. Several limitations of our study should be acknowledged. First, GPs participation was selective. GPs who volunteered to participate in the study could have been especially interested in mental health, however socio-occupational characteristics such as gender, seniority and working time did not differ between GPs answering the Order of Medicine survey and those including patients and responding to the study questionnaire. The GP response rate could have been increased by repeated reminders. The requirement to include only patients consulting for mental distress at work for the first time limited the number of inclusions, but reduced selection bias. The use of the same questions in the Order of Medicine survey and in the study questionnaire allowed sample representativeness and the rate of GP consultation for mental distress at work to be assessed. A patient selection should also be considered, however the physicians participating in the study were instructed to ask all eligible patients to participate in the study. A third limitation was to do with the definitions of work-related psychiatric disorders: the attribution of work-relatedness makes by the GPs or the patient is subjective. For the measurement of GP-diagnosed psychiatric disorders, there were no standardized procedures for the diagnosis. This could have introduced bias in the measurement however in France GPs do not use standardized procedures to establish diag-noses of psychiatric disorders among their patient. Besides we have to cautious about the results because 406 of the small sample size which does not allow the 407 results to be generalized to the French population.

409

  In the present study, most GPs in the Loire area 410 reported difficulties in dealing with distress at work. 411 Screening for the impaired physical and/or psycho-412 logical health status in case of distress at work 413 reported by the patient should guide the GP's clin-414 ical approach. It is not the GP's role to establish a 415 link between symptoms and work. The slowness of 416 reaching a diagnosis of mental distress at work is a 417 real problem: it is alarming to find that patients wait 418 several weeks or even years before consulting their 419 GP, despite experiencing impact on family life and 420 health. For preventive purposes, the factors underly-421 ing such complaints need to be clarified. Referral to 422 occupational medicine could be improved by improv-423 ing communication between patient and GP. It should 424 be borne in mind that all employees have the right to 425 ask for a consultation with their occupational physi-426 cian [17]. The present study highlights the importance 427 of reporting negative experience of working condi-428 tions and impaired mental health to the occupational 429 physician, so that he or she can undertake preventive 430 measures.

were rendered anonymous. Descriptive anal- 124 ysis used SAS software, version 9.3 for chi² or Fisher 125 tests to compare frequencies, with the significance 126 threshold set at p < 0.05. 127 3. Results 128

		99
	-socio-occupational and demographic items: age,	100
	gender, marital status, socio-occupational cate-	101
	gory;	102
	-characteristics of the mental distress: symptom	103
	onset, trigger factors, impact on personal life and	104
	health, level of stress;	105
	-expectations regarding the consultation.	106
	The GP questionnaire covered:	107
	-socio-occupational items: gender, seniority,	108
	number of days worked per week; these data	109
	were collected using the same questions as in the	110
	OM survey of all Loire area GPs, so as to deter-	111
	mine the rate of consultation for mental distress	112
	at work per GP;
		119
	on November 21, 2016, and by University Hospi-	120
	tal ethics committee and the regional review board	121
	(respectively, January 12, 2017, Ref.: IORG0007394,	122
	and June 29, 2017, Ref.: 2017-A01404-49).
	Eight hundred GPs were asked by the local Order	129
	of Medicine to fill out an on-line questionnaire. 118	130
	(71 female, 47 male) did so: i.e., response rate, 14.7%.	131
	Sixteen of the 22 GPs eligible to include patients for	132
	the second study did so, and included 22 patients (15	133
	female, 7 male). Socio-occupational characteristics	134
	such as gender, seniority and working time did not	135
	differ between GPs answering the Order of Medicine	136
	survey and those including patients and responding	137
	to the study questionnaire (Table 1). The mean rate	138
	of consultation for mental distress at work was 2.9%	139
	on the OM survey and 2.1% in the study series.

113

-medical treatment: referral, medical prescription 114 if any, sick-leave if any, diagnosis, etc. 115 The GP questionnaire was completed at the end 116 of the consultation. The study was approved by 117 the French data protection commission (Commission 118 Nationale de l'Informatique et des Libertés: CNIL) 123 Data

Table 1

 1

	Socio-occupational comparison between General Practitioners (GP) sample and Order of
	Medicine (OM) survey of GPs in the Loire area
	Variables		GPs	OM	P-value
			n	%	n	%
	Gender	Male	8	50	47	39.8	0.60
	Female	8	50	71	60.2
	Seniority (years)	<5	2	12.5	30	25.4	0.49
	[5-15]	7	43.7	36	30.5
	[15-25]	4	25.0	21	17.8
		≥ 25	3	18.8	31	26.3
	Half-days worked per week	<7	6	37.5	25	21.2	0.16
		[7-8]	2	12.5	16	13.6
		[8-9]	5	31.3	25	21.2
		≥9	3	18.7	52	44.0
	Patients per half-day	<10	2	12.5	8	6.8	0.36
	[10-14]	7	43.7	58	49.1
	[14-15]	2	12.5	6	5.1
		≥15	5	31.3	46	39.0
	Distress patients per week in last month	0	3	18.7	1	0.8	0.025
		[1-3]	10	62.5	63	53.4
		[3-4]	2	12.5	17	14.4
		≥4	1	6.3	37	31.4
	GPs: General Practitioners sample. OM: Loire Order of Medicine survey.
		Table 2
	Characteristics of General Practitioners' medical treatment
	Variables				n	%
	Longer consultation		No		4	18.2
			Yes		18	81.8
	Psychology referral		No		14	63.6
			Yes		8	36.4
	Psychiatry referral		No		20	90.9
			Yes		2	9.1
	Occupational medicine referral		No		16	72.7
			Yes		6	27.3
	Specialist referral		No		22	100
			Yes		0	0
	Medical treatment		No		12	54.5
			Yes		10	45.5
	Sick-leave certificate		No		0	0
			Yes		22	100
	Work accident certificate		No		22	100
			Yes		0	0
	Occupational disease certificate		No		22	100
			Yes		0	0
	Sick-leave duration		<7 days		1	4.5
			[7-14 days]	11	50.0
			≥ 15 days		10	45.5
	Sick-leave reason		Anxiety		5	22.7
	Anxiety-depressive syndrome	8	36.4
			Depression		1	4.5
			Burn-out		5	22.7
			Asthenia		2	9.1
			Bursitis		1	4.5
	"Reaction" status declared		No		16	72.7
	on medical certificate		Yes		6	27.3
	Medical management of patients consulting for

141 the first time for mental distress at work was char-142 acterized by longer than usual consultation time, 143 systematic sick leave (of less than two weeks in half 144 of the cases), and medical treatment in half of the cases. Reasons for sick leave were mental pathology (anxiety, depression, anxiety-depressive syndrome) in two-thirds of cases and burnout in one-fifth. No

Table 3

 3 Socio-occupational characteristics of patients and GP consultation factors

	Variables

150

referred to an occupational physician (Table

2

). 163 a sick-leave prescription, but less than a quarter 164 expected medical treatment. A minority expected 165 referral to the occupational physician, a psychologist 166

accepted in 596 cases, 40% up on 2015. According to

April and August 2014, gave one of the first estimates

Rehabilitation. 2012;1(22):88-96.

Conflict of interest

The authors declare that they have no competing