Hassan S Bakouch
email: hnbakouch@yahoo.com

Nauman Muhammad

Tassaddaq Khan

Christophe Hussain

A Chesneau

Power

Muhammad Nauman Khan

Tassaddaq Hussain

Christophe Chesneau
email: christophe.chesneau@unicaen.fr

Muhammad Nauman

A Power Log-Dagum Distribution: Estimation and Applications

Keywords: Distributions on the real line, Moments, Estimation, Goodness of fit statistics, TTT-plot. 2000 MSC: 60E05, 62E15

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Statistical distributions play a significant role in describing and predicting real world phenomena. In the 1970s, Camilo Dagum developed a statistical distribution to fit empirical income and wealth data that are not satisfied with the classical distributions (Pareto and lognormal distributions). He looked for a model accommodating the heavy tails appear in empirical income and wealth data distributions, where the former distribution is well captured by the Pareto but not by the lognormal and the latter by the lognormal but not the Pareto. Experimenting with a shifted log-logistic distribution [START_REF] Dagum | A model of income distribution and the conditions of existence of moments of finite order[END_REF], Dagum realized that a further parameter was needed to such distribution which led to the Dagum type I and generalizations with three-parameter and four-parameter distributions [START_REF] Dagum | A new model of personal income distribution: Specification and estimation[END_REF][START_REF] Dagum | The generation and distribution of income, the Lorenz curve and the Gini ratio[END_REF].

In the same era Mielke and Johnson [START_REF] Mielke | Some generalized beta distributions of the second kind having desirable application features in hydrology and meteorology[END_REF] proposed the generalized beta distribution of the second kind abbreviated as GBDII. This distribution is used in the flood fre-quency analysis and it has the beta-k distribution as a sub-model. After that various authors have shown that the Dagum distribution [START_REF] Dagum | A model of income distribution and the conditions of existence of moments of finite order[END_REF] and GBDII are identical and they are two different parameterizations of the same distribution (see, for example, [START_REF] Chotikapanich | Modeling income distributions and Lorenz curves[END_REF][START_REF] Kleiber | Statistical size distribution in economics and actuarial sciences[END_REF]). Domma and Perri [START_REF] Domma | Some developments on the log-Dagum distribution[END_REF] proposed the log-Dagum (LD) distribution obtained by a logarithmic transformation of the Dagum distribution. The LD distribution is defined on the real line and its shape is leptokurtic. Also, it may be symmetric and asymmetric, and hence shall be useful in modeling skewed and leptokurtic distributions which frequently occur in several areas such as finance, reliability, econometrics, insurance and hydrology.

Interpretations of the real world phenomena needs introducing new statistical distributions, namely ones defined on the whole real line and having bimodal behavior for both density and hazard rate functions. Therefore, we introduce a new model called a power log-Dagum (PLD) distribution with the cumulative distribution function (cdf)

F (x) = 1 + e -(ϑx+sign(x) ϑ |x| ϑ ) -ζ , x, ϑ ∈ , ζ > 0, ≥ 0, ( 1
)
where sign is 1 if x > 0, 0 if x = 0 and -1 if x < 0. In addition to, it has the following representation: F (x) = G [x w(x)] where w(x) is the polynomial weight: w(x) = ϑ + ϑ |x| ϑ-1 , satisfying lim x→-∞ xw(x) = -∞ and lim x→+∞ xw(x) = +∞, and G(x) = {1 + e -x } -ζ is a cdf of the LD distribution with parameters (ζ, 1, 1).

The corresponding probability density function (pdf) and the hazard rate function (hrf) are given as

f (x) = ζ ϑ + |x| ϑ-1 e -(ϑx+sign(x) ϑ |x| ϑ ) 1 + e -(ϑx+sign(x) ϑ |x| ϑ ) -(ζ+1) , (2)
h(x) = ζ ϑ + |x| ϑ-1 e -(ϑx+sign(x) ϑ |x| ϑ ) 1 + e -(ϑx+sign(x) ϑ |x| ϑ ) -(ζ+1) 1 -1 + e -(ϑx+sign(x) ϑ |x| ϑ ) -ζ , (3)
respectively.

Obviously, the PLD distribution is defined on the entire real line and this is one of the important features of it, unlike the Dagum [START_REF] Dagum | A model of income distribution and the conditions of existence of moments of finite order[END_REF] and GBDII [START_REF] Mielke | Some generalized beta distributions of the second kind having desirable application features in hydrology and meteorology[END_REF] distributions, which can only provide support on the positive real line.

The PLD distribution defined by (1) has the following submodels.

(1) When ζ = 1, then (1) reduces to the power logistic (PLo) distribution with the density

f (x) = ϑ + |x| ϑ-1 e -(ϑx+sign(x) ϑ |x| ϑ ) 1 + e -(ϑx+sign(x) ϑ |x| ϑ ) 2 , x, ϑ ∈ , ϑ > 0, ≥ 0. ( 4
)
(2) When ϑ = 2, then (1) reduces to the linear log-Dagum distribution with the density

f (x) = ζ (2 + |x|) 1 + e -(2x+sign(x) 2 x 2 ) -ζ-1
e -(2x+sign(x) 2 x 2 ) .

(3) When ϑ = 2 and ζ = 1, then [START_REF] Aarset | How to identify bathtub hazard rate[END_REF] gives the linear logistic distribution with the density

f (x) = (2 + |x|) e -(2x+sign(x) 2 x 2 ) 1 + e -(2x+sign(x) 2 x 2 )
2 .

(4) When ϑ = 2 and = 0, then (1) reduces to the log-Dagum distribution with parameters (ζ, 2) and the density

f (x) = 2ζ 1 + e -2x -ζ-1 e -2x ,
which is introduced by [START_REF] Domma | Some developments on the log-Dagum distribution[END_REF]. ( 5) When ϑ = 2, ζ = 1 and = 0, then (1) leads to the known logistic distribution with the density

f (x) = 2e -2x [1 + e -2x ] 2 .
[Figure 1 [Figure 2 about here.]

Plots of Figure 2 display the density functions of the PLD distribution. Figure 2 portrays that changing against the fixed ϑ and ζ shift the mode towards left. But in case of changing ζ with fixed ϑ and shift the curve towards right. However, the interesting feature of the distribution is its bimodal behavior which is frequently used in biomedical and engineering phenomenon like formation of bathtub shapes of the hazard function. Such bimodal behavior is captured while fixing ζ and and changing ϑ. Further, Figure 2 portrays that bimodality nature of the curve shift towards right as ϑ increases.

[Figure 3 about here.]

Hazard function is an important indicator for observing the deteriorating condition of a product which ranges from increasing, decreasing, bathtub (BT) to inverse bathtub (IBT) shapes. So in this regard Figure 3 speaks out it self and justifies the potential of the model. Moreover, the hazard function plots in Figure 3 also portray the deteriorating conditions of the product as time increases in terms of spontaneous spikes at the end of either increasing or decreasing hazard rate. This implies that the hazard function is sensitive against different combinations of the parameters as time changes, which seems to be a refine image of non stationarity process and hence the hazard curve does not remain stable as times passes. Moreover, Figure 3 displays increasing, decreasing, bimodal and upside down bathtub hazard shapes.

It is worth mentioning that the economic and hydrologic data analysis is based on the assumption that the data are stationary but a number of documented researches showed that these data may be non-stationary. Therefore, identification and use of non-stationary probabilistic models in practice have been recommended (See [START_REF] Salas | Revisiting the concepts of return period and risk for non-stationary hydrologic extreme events[END_REF]). In application section of this paper, Example 2 shows the merit of the PLD model for modeling and analyzing non-stationary time series data as inflation rates with positive and negative values.

The reminder of the paper is outlined as follows. We discuss some statistical properties of the PLD distribution in Section 2, including moments, moment generating function and moments of order statistics. We provide in Section 3 three estimation procedures of the PLD model parameters, namely the maximum likelihood estimation, ordinary and weighted least square estimation, and they are compared using simulations studies. Applications of the PLD model for three practical data sets are discussed in Section 4.

Some Statistical Properties

In this section, we study some statistical properties of the PLD distribution, including moments, moment generating function and moments of order statistics.

Moments and moment generating function

Let X be a PLD random variable, the r th moment of X, say µ r = E(X r ), follows from (2) as

µ r = ζ ∞ -∞ x r ϑ + |x| ϑ-1 e -(ϑx+sign(x) ϑ |x| ϑ ) 1 + e -(ϑx+sign(x) ϑ |x| ϑ ) -ζ-1 dx.
Using the series expansion

1 + e -(ϑx+sign(x) ϑ |x| ϑ ) -ζ-1 =              ∞ j=0 -ζ -1 j e -j(ϑx+ ϑ x ϑ ) , if x > 0, ∞ j=0 -ζ -1 j e (ζ+j+1)(ϑx-ϑ (-x) ϑ ) , if x < 0, (5)
we obtain the moment as

µ r = ζ ϑ ∞ j=0 -ζ -1 j 0 -∞ x r e -(ϑx-ϑ (-x) ϑ ) e (ϑx-ϑ (-x) ϑ )(ζ+j+1) dx + ∞ 0 x r e -(ϑx+ ϑ x ϑ ) e -j(ϑx+ ϑ x ϑ ) dx + ζ ∞ j=0 -ζ -1 j 0 -∞ x r (-x) ϑ-1 e -(ϑx-ϑ (-x) ϑ ) e (ϑx-ϑ (-x) ϑ )(ζ+j+1) dx + ∞ 0
x r+ϑ-1 e -(ϑx+ ϑ x ϑ ) e -j(ϑx+ ϑ x ϑ ) dx .

By using the series expansion, we can write

µ r = ζ ϑ ∞ j,k=0 -ζ -1 j (-ϑ (ζ + j)) k k! 0 -∞ x r (-x) ϑ k e ϑ(ζ+j)x dx + (-ϑ (1 + j)) k k! ∞ 0 x r+ϑ k e -ϑ(1+j)x dx + ζ ∞ j,k=0 -ζ -1 j (-ϑ (ζ + j)) k k! 0 -∞ x r (-x) ϑ(1+k) -1 e ϑ(ζ+j)x dx + (-ϑ (1 + j)) k k! ∞ 0
x r+ϑ(1+k) -1 e -ϑ(1+j)x dx .

Hence

µ r = ζ ϑ ∞ j,k=0 -ζ -1 j 1 k! {-ϑ (ζ + j)} k (-1) r {ϑ(ζ + j)} 1+ϑk+r + {-ϑ (1 + j)} k (ϑ + ϑj) 1+ϑk+r Γ (1 + ϑk + r) + ζ ∞ j,k=0 -ζ -1 j 1 k! {-ϑ (ζ + j)} k (-1) r {ϑ(ζ + j)} ϑ+ϑk+r + {-ϑ (1 + j)} k (ϑ + ϑj) ϑ+ϑk+r Γ (ϑ + ϑk + r) .
The moment generating function (mgf) of the PLD distribution is

M (t) = ζ ∞ -∞ e tx ϑ + |x| ϑ-1 1 + e -(ϑx+sign(x) ϑ |x| ϑ ) - (ζ+1)
e -(ϑx+sign(x) ϑ |x| ϑ ) dx, (6) and making use of (5), we have

M (t) = ζ ϑ ∞ j,k=0 -ζ -1 j 1 k! {-ϑ (ζ + j)} k {ϑ(ζ + j) + t} 1+ϑk + {-ϑ (1 + j)} k (ϑ + ϑj -t) 1+ϑk Γ (1 + ϑk) + ζ ∞ j,k=0 -ζ -1 j 1 k! {-ϑ (ζ + j)} k {ϑ(ζ + j) + t} ϑ(1+k) + {-ϑ (1 + j)} k (ϑ + ϑj -t) ϑ(1+k) Γ (ϑ + ϑk) .

Moments of order statistics

Let X 1 , X 2 , . . . , X n be a random sample from the PLD distribution and its order statistics is X 1:n , X 2:n , . . . , X n:n . Let f i:n (x) and F i:n (x) denote the pdf and the cdf of the i th order statistic X i:n , respectively. Hence, using the standard expressions of order statistics, we find that

f i:n (x) = n! (i -1)! (n -i)! n-i l=0 n -i l (-1) l [F (x)] i-1+l f (x) ,
and

F i:n (x) = n! (i -1)! (n -i)! n-i l=0 n -i l (-1) l i + l [F (x)] i+l ,
where F (x) and f (x) are given by equations ( 1) and ( 2), respectively. Hence, the r th moment of the i th order statistic X i:n is

E (X r i:n ) = ∞ -∞ x r f i:n (x) dx,
and by making use of f i:n (x) and equation ( 5) with some algebra, one find that

E (X r i:n ) = ζ n! (i -1)! (n -i)! n-i l=0 ∞ j,k=0 n -i l -ζ(i + l) -1 j (-1) l k! × ϑ { ϑ (1 + ζ(i + l) -j)} k (-1) r {ϑ(j -1 -ζ(i + l))} 1+ϑk+r + {-ϑ (1 + j)} k (ϑ + ϑj) 1+ϑk+r Γ (1 + ϑk + r) + { ϑ (1 + ζ(i + l) -j)} k (-1) r {ϑ(j -1 -ζ(i + l))} ϑ+ϑk+r + {-ϑ (1 + j)} k (ϑ + ϑj) ϑ+ϑk+r Γ (ϑ + ϑk + r) .

Estimation Procedures

In this section we investigate three estimation procedures for the PLD model parameters, namely the maximum likelihood estimation, ordinary and weighted least square estimations. Performance of those procedures is investigated by simulation studies.

Maximum likelihood estimators

Let X 1 , X 2 , • • • , X n be a random sample from the PLD distribution with parameter vector Θ=(ϑ, ζ, ) and x 1 , x 2 , • • • , x n are the corresponding observed values, then the joint probability function of X 1 , X 2 , • • • , X n as a log-likelihood function can be expressed as

(Θ) = n log (ζ) + n i=1 log ϑ + |x i | ϑ-1 - n i=1 ϑx i + sign(x i ) ϑ |x i | ϑ -(ζ + 1) n i=1 log 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ . (7)
The associated nonlinear log-likelihood equations ∂ (Θ) ∂Θ = 0 are given by

∂ (Θ) ∂ϑ = n i=1 1 + log (|x i |) |x i | -1+ϑ ϑ + |x i | -1+ϑ - n i=1 x i + ϑ sign(x i )|x i | ϑ log(|x i |) - 1 ϑ -(1 + ζ) n i=1 e -ϑxi-sign(xi) ϑ |xi| ϑ -x i -ϑ sign(x i )|x i | ϑ log(|x i |) -1 ϑ 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ = 0, ∂ (Θ) ∂ζ = n ζ - n i=1 log 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ = 0, ∂ (Θ) ∂ = - n i=1 sign(x i ) |x i | ϑ ϑ + (1 + ζ) n i=1 e -ϑxi-sign(xi) ϑ |xi| ϑ sign(x i )|x i | ϑ 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ ϑ + n i=1 |x i | -1+ϑ ϑ + |x i | -1+ϑ = 0.
Solving the equations above simultaneously, we obtain the maximum likelihood estimators (MLEs) of the model parameters.

The information matrix is used to establish the confidence intervals for the model parameters. The elements of the 3 × 3 observed information matrix J(Θ) = J rs (Θ) for r, s = ϑ, ζ, are provided in Appendix B.

In the simulation study for the MLEs (Subsection 3.3), it is observed that the bias and mean square error of the MLEs decrease as sample sizes increase which suit the usual norms of asymptotic properties of the MLEs. However, when all the parameters posses values greater than one, a downward bias is only reflected with ϑ and , and when all the parameters posses values less than one, a downward bias is only reflected with ϑ. In addition to, an upward bias is observed for all other combinations. Those features can be observed in Table 1.

Ordinary and weighted least-square estimators

Let x 1 , x 2 , • • • , x n be an ordered sample of the random sample of size n from a distribution function F (•). Let F (•) follow the PLD distribution, then the least square estimators (LSE) can be obtained by minimizing

L(Θ) = n i=1 1 + e -ϑ xi-sign(xi) ϑ |xi| ϑ -ζ - i n + 1 2 , (8)
with respect to the unknown parameters of the distribution. The associated nonlinear equations ∂L(Θ) ∂Θ = 0 are given by

∂L(Θ) ∂ϑ = 2ζ n i=1 e -ϑxi-sign(xi) ϑ |xi| ϑ 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -ζ - i 1 + n × 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -1-ζ x i + ϑ sign(x i )|x i | ϑ log(|x i |) - 1 ϑ = 0, ∂L(Θ) ∂ζ = -2 n i=1 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -ζ 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -ζ - i 1 + n × log 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ = 0, ∂L(Θ) ∂ = 2ζ ϑ n i=1 e -ϑxi-sign(xi) ϑ |xi| ϑ 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -1-ζ × 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -ζ - i 1 + n sign(x i )|x i | ϑ = 0.
Solving this system of nonlinear equations simultaneously will yield the LSEs of the distribution parameters.

In the simulation study, a down bias is noted in all of the parameter combination except for , when it is less than one. Moreover, fluctuation in bias and mean square error is reflected in parameter combinations which obviously approaches to zero. Table 2 shows such reflections.

On the other side, the weighted least square estimators (WLSE) of parameters of the PLD distribution can be obtained by minimizing

W(Θ) = n i=1 (n + 1) 2 (n + 2) i(n -i + 1) 1 + e -ϑ xi-sign(xi) ϑ |xi| ϑ -ζ - i n + 1 2 , (9)
with respect to the unknown parameters. The associated nonlinear equations ∂W(Θ) ∂Θ = 0 are given by

∂W(Θ) ∂ϑ = 2ζ n i=1 (n + 1) 2 (n + 2) i(n -i + 1) e -ϑxi-sign(xi) ϑ |xi| ϑ 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -1-ζ × 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -ζ - i 1 + n × x i + ϑ sign(x i )|x i | ϑ log(|x i |) - 1 ϑ = 0, ∂W(Θ) ∂ζ = -2 n i=1 (n + 1) 2 (n + 2) i(n -i + 1) 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -ζ - i 1 + n × 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -ζ log 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ = 0, ∂W(Θ) ∂ = 2ζ ϑ n i=1 (n + 1) 2 (n + 2) i(n -i + 1) e -ϑxi-sign(xi) ϑ |xi| ϑ 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -1-ζ × 1 + e -ϑxi-sign(xi) ϑ |xi| ϑ -ζ - i 1 + n sign(x i )|x i | ϑ = 0.
Solving the system of equations above yields the WLSEs of the distribution parameters.

A simulation work, see the next subsection, shows that a downward bias is observed in all parameter combinations except when < 1. In addition, fluctuation in bias and mean square error is remarked in all parameter combinations which approaches to zero. Table 3 exhibits such reflections.

Simulation studies

In this subsection, we perform a small simulation study to gain some results on the estimators obtained by using the above estimation methods as well as their asymptotic behavior for finite samples. All the numerical computations are performed via Mathematica 8.0 using the random numbers generator code. We consider the following different • Average mean square error (MSE) of the simulated estimates:

1 1000 1000 i=1 (Θ -Θ) 2 ,
where Θ = ( θ, ˆ , ζ) are estimates of the parameter vector Θ = (ϑ, ζ, ). The results are reported in Tables 1, 2 and 3.

[Table 1

Evaluation Measures and Practical Data Examples

In the statistical literature, various distributions are employed for different situations. Some of them are suitable for symmetrical data, some for skewed data and others are for both aspects. In this regard, we study those distributions which are either, positively skewed, negatively skewed or symmetrical and posses increasing, decreasing, bathtub shapes (BTS) and inverse bathtub (IBT) hazard function. Moreover, comparison also requires that random variable should also be defined on the whole real line, i.e. x ∈ . The above conditions demand the following known distributions that are: Normal, Cauchy, Gumbel, Dagum. For the probability density functions of the above mentioned distributions, the reader may refer to [START_REF] Johnson | Continuous univariate distributions Volume-1 Second edition[END_REF] and [START_REF] Johnson | Continuous univariate distributions Volume-2 Second edition[END_REF], noting that the log-dagum (LD) was developed by [START_REF] Domma | Some developments on the log-Dagum distribution[END_REF].

Numerical measures

In order to demonstrate the proposed methodology, we consider three different practical data sets described below with their analysis. They represent different level of skewness ranging from negative skewness to positive skewness along with various demonstration of failure rate pattern, like increasing, decreasing, and bathtub shape.

Moreover, perfection of competing models is also tested via the chi-squared (χ 2 ), Kolmogrov-Simnorov(K S), the Anderson-Darling (A * ) and the Cramér-von Misses (W * ) statistics. The mathematical expressions for the statistics are given by

KS = max i m -z i , z i - i -1 m , χ 2 = m i=1 (o i -e i ) 2 e i , A * = 2.25 m 2 + 0.75 m + 1 -m - 1 m m i=1 (2i -1) ln(z i (1 -z m-i+1 )) , W * = m i=1 z i - 2i -1 2m 2 + 1 12m
,

where m denotes the number of classes, z i = F X (x i ), the x i 's being the ordered observations, o i and e i are the observed and expected frequencies of the i th class, respectively.

Graphical measure

A hazard or failure rate function can be considered for finding the chances of occurring a dangerous event that can lead to an emergency or disaster. It also conveys the meaning of intensity function, or risk rate, among other names. Moreover, such function also portrays different meaning of inherent characteristics of different lifetime phenomena, like increasing hazard rate (IHR), constant, decreasing hazard rate (DHR), bathtub (BT) and inverse bathtub (IBT) shapes. For modeling the above mentioned characteristic of the hazard function, researchers generally used a graphical tool known as the total time on test (TTT) plot. The TTT-plot is a graph which mainly serves to discriminate between different types of aging represented in hazard rate shapes, for details the readers are referred to [START_REF] Rinne | The Weibull distribution; A Handbook[END_REF]. The TTT plot is drawn by plotting T ( i n ) = i r=1 y r:n + (n -i)y i:n n r=1 y r:n against i n , where i = 1, . . . , n and y r:n , r = 1, . . . , n are the order statistics of the sample. An empirical TTT-plot that takes its course randomly around the 45 o -line indicates a sample from an exponential distribution, TTT-plot that is nearly concave (convex) and is mainly above (below) the 45 o -line indicates a sample from an IHR (DHR) distribution, and TTT-plot that takes concave (convex) and then convex (concave) is related to a BT (IBT) hazard rate. The TTT-plots for the next data sets are displayed in Figure 4, which reveal increasing IHRs for the first two data sets while the third data set depicts the BT hazard rate.

[Figure 4

about here.]

Example 1. The first data set represents the strength data originally reported in [START_REF] Badar | Statistical aspects of fiber and bundle strength in hybrid composites[END_REF]. It represents the strength measured in GPA for single carbon fibers and impregnated 1000-carbon fiber tows. Carbon fiber tow is a bunch of carbon fibre, like a yarn, that is commonly woven into carbon fibre fabrics. Single fibers were tested under tension at gauge length of 10 mm. Such data have been analyzed previously by [START_REF] Gupta | Generalized logistic distributions[END_REF] and [START_REF] Asgharzadeh | A generalized skew logistic distribution[END_REF]. The data are as follows: 3.145, 3.220, 3.223, 3.235, 3.243, 3.264, 3.272, 3.294, 3.332, 3.346, 3.377, 3.408, 3.435, 3.493, 3.501, 3.537, 3.554, 3.562, 3.628, 3.852, 3.871, 3.886, 3.971, 4.024, 4.027, 4.225, 4.395, 5.020.

We fitted all six distributions to the above data by the method of maximum likelihood. Table 4 presents the descriptive statistics. Table 5 portrays parameter estimates, χ 2 , KS with p -values Anderson-Darling and Cramér-von Misses statistics.

[Table 4 about here.]

[Table 5

about here.]

Analysis: From Tables 4, and 5, we see that the proposed PLD model is the most suitable model when the data is positively skewed and underdispersed with approximately equal mean deviation about mean and median as well as posses IFR behavior. Moreover, the goodness of fit statistics also show that the proposed model is the most promising one in such situations too. In addition to, the proposed model also shows promising results while comparing with the models discussed by [START_REF] Asgharzadeh | A generalized skew logistic distribution[END_REF] in terms of the smallest χ 2 value. The compared models are: logistic distribution L(µ, σ) with χ 2 = 0.8832, the type III generalized logistic distribution GL(µ, σ, α) with χ 2 = 1.1055, the skew logistic distribution SL(µ, σ, λ) with χ 2 = 0.8345, generalized skew logistic GSL(µ, σ, λ, α) with χ 2 = 0.2682 and Azzalini and Capitanio's skew t distribution ST(µ, σ, λ, α) with χ 2 = 0.2684. For the density functions of the discussed distributions the reader may refer to [START_REF] Asgharzadeh | A generalized skew logistic distribution[END_REF] and the references therein. Moreover, while comparing with [START_REF] Gupta | Generalized logistic distributions[END_REF] proportional reversed hazard logistic (PRHL) and skew logistic distributions, it is observed that PLD has the smallest χ 2 and KS values.

Example 2. The second data set displays the inflation rate of USA for the period 1965 to 1981, and it is defined as the rate at which prices increase over time, resulting in a fall in the purchasing value of currency. The data is reported by [START_REF] Gujarati | Basic econometrics, Fourth Edition[END_REF], pp. 502. The data values are: -0.4, 0.4, 2.9, 3.0, 1.7, 1.5, 1.8, 0.8, 1.8, 1.6, 1.0, 2.3, 3.2, 2.7, 4.3, 5.0, 4.4, 3.8, 3.6, 7.9, 10.8, 6.0, 4.7, 5.9, 7.9, 9.8, 10.2, 7.3.

We fitted all of the distributions except Dagum to the above data by the method of maximum likelihood. Table 6 identifies the descriptive statistics behavior, Table 7 presents the parameter estimates, χ 2 , p -values, KS, A * and W * .

[Table 6 about here.] [Table 7 about here.] Analysis: From Tables 6 and7, we observe that, the proposed model is also suitable when data are overdispersed and based on negative values. Also, the model has positive skewness and demonstrates IFR failure rates which can be visualized in TTT plot of data II. In addition to, the proposed model is also one of the appropriate model from the perspective of goodness of fit statistics. The importance of model is also identified in least χ 2 and the largest p -value.

[Figure 5 about here.] From econometrics point of view, the sample path of the data given by Figure 5 shows that the data series are non-stationary. Moreover, the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test rejects the stationarity (p-value < 0.05) of the data. Also, following the Box-Ljung test, there is a structure of dependence behind those time series data (p -value < 0.001).

[Figure 6 about here.] Example 3. The third data set is the famous Aarset data which gives the uncensored practical data set and it refers to the lifetimes of 50 components reported by [START_REF] Aarset | How to identify bathtub hazard rate[END_REF]. The values of the data are: 0.1, 0.2, 1, 1, 1, 1, 1, 2, 3, 6, 7, 11, 12, 16, 18, 18, 18, 18, 18, 20, 21, 32, 36, 40, 46, 47, 50, 55, 60, 63, 63, 67, 67, 67, 67, 72, 75, 79, 82, 83, 84, 84, 84, 85, 85, 85, 85, 85, 86, 86. We fitted all six distributions to the data above by the method of maximum likelihood. Table 8 gives the descriptive statistics of them, Table 9 presents the parameter estimates, χ 2 , p -values, KS, A * and W * .

[Table 8 about here.]

[Table 9 about here.] Analysis: From Table 8, the data are overdispersed and negatively skewed. Moreover, the TTT plot for the Aarset's data indicates BT shapes. In Table 9, we have observed that PLD is the most suitable model based on the goodness of fit statistics with p-vaalue. Further, Figure 6 also demonstrates the flexibility of PLD distribution.

Appendices

Appendix A. Algorithm: i Generate U i , i = 1, 2, 3, ..., n from uniform distribution on the interval (0,1).

ii Generate ln U

-1 ζ i -1 -ϑ , i = 1, 2, 3, ..., n.
iii Set

E i = ln U -1 ζ i -1 -ϑ ϑ+ 2 , i = 1, 2, 3, ..., n. iv Set D i = ln U -1 ζ i -1 -ϑ 1 ϑ ϑ+ 2 , i = 1, 2, 3, ..., n. v Set Y 1 = E i , i = 1, 2, 3, ..., n. vi Set Y 2 = D i , i = 1, 2, 3, ..., n. ⇒ If U i ≤ ϑ 2 + ϑ 2 , then set X i = Y 1 otherwise set X i = Y 2 , i = 1, 2, 3, ..., n. J ϑ ϑ (Θ) = n i=1    |x i | ϑ-1 log (|x i |) 2 ϑ + |x i | ϑ-1 - 1 + |x i | ϑ-1 log (|x i |) 2 ϑ + |x i | ϑ-1 2    - n i=1 2 -2 ϑ log (|x i |) + ϑ 2 log (|x i |) 2 sign (x i ) |x i | ϑ ϑ 3 -(1 + ζ) n i=1 -e -2ϑxi-2 sign(x i) |x i | ϑ ϑ 1 + e ϑxi+ sign(x i) |x i | ϑ ϑ ϑ × 2 -2 ϑ log (|x i |) + ϑ 2 log (|x i |) 2 sign (x i ) |x i | ϑ + ϑ 2 x i -sign (x i ) |x i | ϑ (1 -ϑ log (|x i |)) 2 -1 + e ϑxi+ sign(x i) |x i | ϑ ϑ × ϑ 2 x i -sign (x i ) |x i | ϑ (1 -ϑ log (|x i |)) 2 × 1 1 + e -ϑxi-sign(x i) |x i | ϑ ϑ 2 ϑ 4 , J ϑ ζ (Θ) = n i=1 ϑ 2 x i -sign (x i ) |x i | ϑ + ϑ log (|x i |) sign (x i ) |x i | ϑ 1 + e ϑxi+ sign(x i) |x i | ϑ ϑ ϑ 2 , J ζ ζ (Θ) = - n ζ 2 , J ζ (Θ) = n i=1 - e -|xi| ϑ sign(x i) ϑ -ϑ xi |x i | ϑ sign (x i ) ϑ 1 + e -|xi| ϑ sign(x i) ϑ -ϑ xi , J ϑ (Θ) = - n i=1 (-1 + ϑ log (|x i |)) sign (x i ) |x i | ϑ ϑ 2 + n i=1 (-1 + ϑ log (|x i |)) sign (x i ) x i |x i | ϑ (ϑx i + sign (x i ) |x i | ϑ ) 2 -(1 + ζ) n i=1 sign (x i ) |x i | ϑ ϑ + e ϑxi+ sign(x i) |x i | ϑ ϑ ϑ -ϑ 2 log (|x i |) -e ϑxi+ sign(x i) |x i | ϑ ϑ ϑ 2 log (|x i |) +e ϑxi+ sign(x i) |x i | ϑ ϑ ϑ 2 x i -e ϑxi+ sign(x i) |x i | ϑ ϑ sign (x i ) |x i | ϑ +e ϑxi+ sign(x i) |x i | ϑ ϑ ϑ log (|x i |) sign (x i ) |x i | ϑ × 1 ϑ 3 1 + e ϑxi+ sign(x i) |x i | ϑ ϑ 2 , J (Θ) = n i=1 - x 2ϑ-2 i ϑ + |x i | ϑ-1 2 -(1 + ζ) n i=1      - e -2 |xi| ϑ sign(x i) ϑ -2ϑxi x 2ϑ i ϑ 2 1 + e -|xi| ϑ sign(x i) ϑ -ϑxi 2 + e -|xi| ϑ sign(x i) ϑ -ϑxi x 2ϑ i ϑ 2 1 + e -|xi| ϑ sign(x i) ϑ -ϑ xi     ,
the remaining elements follow by symmetry.

Figure 1

 1 Figure 1 gives the plots of the cumulative distribution function of the PLD distribution. The plots of this figure shows that for fixed ϑ and and changing ζ the curve stretch out insignificantly towards right as ζ increases. However, for fixed ζ and and changing ϑ the curve stretch out towards right significantly as ϑ increases.

  model parameters: Model-I: ϑ = 0.3642, = 2.3532, ζ = 1.2315, Model-II: ϑ = 5.8759, = 5.4314, ζ = 4.5472, Model-III: ϑ = 1.2536, = 0.4721, ζ = 1.1973 and Model-IV: ϑ = 0.0542, = 0.0135, ζ = 0.0432 for all the above discussed methods. We consider the following sample sizes n = 15, n = 25 (small), 50 (moderate), and 150 (large). For each model parameters and for each sample size, we compute the MLEs, LSEs and WLSEs of each ϑ, , ζ. We repeat this process 1000 times and compute the average bias and mean square error (MSE) for all replications in the relevant sample sizes. The analysis computes the coming values: • Average bias of the simulated estimates:

Figure 2 .Figure 3 .

 23 Figure 2.: Density plots of PLD

  1.901, 2.132, 2.203, 2.228, 2.257, 2.350, 2.361, 2.396, 2.397, 2.445, 2.454, 2.474, 2.518, 2.522, 2.525, 2.532, 2.575, 2.614, 2.616, 2.618, 2.624, 2.659, 2.675, 2.738, 2.740, 2.856, 2.917, 2.928, 2.937, 2.937, 2.977, 2.996, 3.030, 3.125, 3.139,

Table 1 .

 1 : Average Bias and MSE Values of MLEs from Simulation of the PLD Distribution

	Parameter	Sample Size Bias( θ) Bias(ˆ ) Bias( ζ) MSE( θ) MSE(ˆ ) MSE( ζ)
	ϑ = 0.3642 = 2.3532 ζ = 1.2315	n = 15 n = 25 n = 50 n = 150	3.2184 3.0733 3.0236 2.9609	2.6517 2.1018 2.0380 1.4949	5.7929 4.6184 3.2913 3.1971	11.6709 10.2525 10.0582 9.4931	13.4664 6.7234 6.1874 2.4417	49.1715 24.3536 11.189 10.4433
	ϑ = 5.8759 = 5.4314 ζ = 4.5472	n = 15 n = 25 n = 50 n = 150	-3.0702 -3.1834 -3.3497 -3.2656	-3.9940 -3.5111 -3.7332 -3.6214	9.4900 9.3246 3.4478 3.3609	10.0832 10.0536 11.2945 10.7069	16.0354 13.5562 13.9725 13.1413	536.739 474.913 15.2051 13.1649
	ϑ = 1.2536 = 0.4721 ζ = 1.1973	n = 15 n = 25 n = 50 n = 150	0.3297 0.3081 0.1806 0.1768	1.1955 0.8850 0.8443 0.6552	0.2712 0.2221 0.1125 0.1088	0.2526 0.2371 0.1390 0.1577	1.5412 0.9822 0.9056 0.5455	0.2954 0.1724 0.0798 0.0534
	ϑ = 0.0542 = 0.0135 ζ = 0.0432	n = 15 n = 25 n = 50 n = 150	-0.0476 -0.0448 -0.0483 -0.0494	0.0856 0.0750 0.0647 0.0629	0.0509 0.0332 0.0143 0.0114	0.0023 0.0020 0.0023 0.0024	0.0080 0.0058 0.0042 0.0039	0.0027 0.0016 0.0003 0.0001

Table 2 .

 2 : Average Bias and MSE Values of LSEs from Simulation of the PLD Distribution

	Parameter	Sample Size Bias( θ) Bias(ˆ ) Bias( ζ) MSE( θ) MSE(ˆ ) MSE( ζ)
	ϑ = 0.3642 = 2.3532 ζ = 1.2315	n = 15 n = 25 n = 50 n = 150	-0.2599 -0.3037 -0.3037 -0.2993	-0.3483 -0.1527 -0.1527 -0.0525	-0.1233 -0.1233 -0.1233 -0.0206	0.0676 0.0922 0.0922 0.0896	0.1213 0.0233 0.0233 0.0027	0.0152 0.0152 0.0152 0.0004
	ϑ = 5.8759 = 5.4314 ζ = 4.5472	n = 15 n = 25 n = 50 n = 150	-5.7381 -5.7627 -5.7738 -5.7657	-1.1245 -1.1072 -1.1013 -0.9983	-0.3368 -0.3363 -0.3363 -0.2261	32.9289 33.2108 33.3378 33.2433	1.2645 1.2261 1.2129 0.9967	0.1134 0.1131 0.1131 0.0511
	ϑ = 1.2536 = 0.4721 ζ = 1.1973	n = 15 n = 25 n = 50 n = 150	-1.2534 -1.2522 -1.1885 -1.224	-0.3868 -0.0672 0.4087 0.4085	-0.9971 -0.0890 -0.1158 -0.1316	1.5711 1.5682 1.4252 1.4984	0.1496 0.0045 0.1836 0.1832	0.9942 0.0079 0.0168 0.0181
	ϑ = 0.0542 = 0.0135 ζ = 0.0432	n = 15 n = 25 n = 50 n = 150	-0.0517 -0.0541 -0.0541 -0.0541	0.0786 0.0255 0.0210 0.0158	-0.0166 -0.0015 -0.0218 -0.0245	0.0026 0.0029 0.0029 0.0029	0.0213 0.0028 0.0023 0.0020	0.0014 0.0020 0.0016 0.0013

Table 3 .

 3 : Average Bias and MSE Values of WLSEs from Simulation of the PLD Distribution

	Parameter	Sample Size Bias( θ) Bias(ˆ ) Bias( ζ) MSE( θ) MSE(ˆ ) MSE( ζ)
	ϑ = 0.3642 = 2.3532 ζ = 1.2315	n = 15 n = 25 n = 50 n = 150	-0.3104 -0.3076 -0.3140 -0.3027	-0.1396 -0.0296 -0.0311 -0.0008	-0.1294 -0.1294 -0.1293 -0.0212	0.0963 0.0947 0.0986 0.0917	0.0195 0.0009 0.0009 0.0000	0.0167 0.0167 0.0167 0.0005
	ϑ = 5.8759 = 5.4314 ζ = 4.5472	n = 15 n = 25 n = 50 n = 150	-5.7403 -5.7528 -5.7544 -5.7478	-0.0774 -0.0602 -0.0792 -0.0661	-0.3370 -0.3369 -0.3369 -0.1261	32.9523 33.0951 33.1141 33.0376	0.0060 0.0034 0.0062 0.0043	0.1136 0.1135 0.1135 0.0159
	ϑ = 1.2536 = 0.4721 ζ = 1.1973	n = 15 n = 25 n = 50 n = 150	-1.2404 -1.2522 -1.1868 -1.2520	0.3802 -0.0671 0.3061 -0.2662	0.0028 -0.0890 -0.1087 -0.1951	1.5386 1.5682 1.4212 1.5677	0.1446 0.0045 0.1772 0.0714	0.0000 0.0079 0.0191 0.0380
	ϑ = 0.0542 = 0.0135 ζ = 0.0432	n = 15 n = 25 n = 50 n = 150	-0.0512 -0.0541 -0.0541 -0.0541	0.0510 0.0374 0.0174 0.0018	-0.0127 -0.0417 -0.0316 -0.0021	0.0026 0.0029 0.0029 0.0029	0.0097 0.0024 0.0004 0.0003	0.0012 0.0017 0.0019 0.0010

Table 4 .

 4 : Descriptive Statistics for Strength Measured in GPA Data

	Sample Size	Mean	Median Variance Skewness Kurtosis	M D μ	M D μ	Entropy
	63	3.0593	2.996	0.3855	0.6328	3.2863	0.5023 0.5023	4.1211

Table 5 .

 5 : Maximum Likelihood Estimates and Goodness of Fit Statistics for Strength Measured in GPA Data

	Distribution	Estimates			χ 2	p -value	K S	A *	W *
	Dagum(µ, σ, λ)	1.9907	7.20072	2.6203	0.6754	0.879	0.1011	0.1028	0.5378
	Cauchy(µ, σ)	2.9902	0.4032		3.1015	0.376	0.1236	0.1804	1.5100
	Normal(µ, σ)	3.0593	0.6159		1.1535	0.764	0.1046	0.0885	0.5985
	Gumbel(µ, σ)	3.3835	0.6993		7.6350	0.054	0.1376	0.2795	1.9607
	LD(ϑ, , ζ)	2.0640	32.6578	10.0233	0.2820	0.991	0.0948 0.08707 0.4423
	PLD(ϑ, , ζ)	1.5238	0.2525	151.0774 0.2583	0.992	0.0945	0.0822	0.4166

Table 8 .

 8 : Descriptive Statistics for Aarset Data

	Sample Size	Mean	Median Variance Skewness Kurtosis	M D μ	M D μ	Entropy
	50	43.866	46.5	1079.3	-0.0224	1.3776	30.0594 30.0594	3.1969

Appendix B. Information Matrix

The 3 × 3 total information matrix along with its elements are given by


Pierre Colin
email: pierre.colin@sano.com

Sandrine Micallef

Maud Delattre

Éric Pierre Mancini

P Colin Ab

M Delattre Bc

E Parent

Towards using a full spectrum of early clinical trial data: a retrospective analysis

Keywords: Late Toxicities, Longitudinal data, Dose-Finding, Oncology

Introduction & Motivation

The nal goal of early phase oncology trials consists in choosing a Recommended Phase 2 Dose (RP2D) among a set of dierent explored doses. Selecting the RP2D is of major interest since it conditions the following clinical development phases [START_REF] Booth | On behalf of the Task Force on Methodology for the Development of Innovative Cancer Therapies, Endpoints and other considerations in phase I studies of targeted anticancer therapy: Recommendations from the task force on Methodology for the Development of Innovative Cancer Therapies (MDICT)[END_REF][START_REF] Ratain | Recommended changes to oncology clinical trial design: revolution or evolution?[END_REF][START_REF] Arrowsmith | Phase II failures: 20082010[END_REF]. On top of the importance of the task, recommending a dose is dicult for many reasons: limited amount of data, heterogeneous population, complex mechanism of distribution and target reaching. Because of study duration constraints or incompleteness of the data, designs generally rely on early acute toxicity but nal dose recommendation should handle any kind of clinical meaningful toxic event and not just the early acute adverse events.

In this paper, we propose an approach reassessing the Maximum Tolerated Dose (MTD) at the end of the trial when most of the late toxicity data are observed. We think that the nal recommended dose should not necessarily be the latest explored one, as it is the case with classical design. Indeed, the latest explored dose usually corresponds to the one closest of a targeted rate of early severe toxicity [START_REF] O'quigley | Continual reassessment method: a practical design for phase I clinical trials in cancer[END_REF], or to the highest limiting early acute toxicity risks [START_REF] Babb | Cancer Phase I clinical trials: Ecient dose escalation with overdose control[END_REF][START_REF] Neuenschwander | Critical aspects of the Bayesian approach to phase I cancer trials[END_REF] but is not necessarily admissible in terms of overall safety.

For more than a decade, targeted therapies have been appearing with complex mechanisms of action targeting specic pathways, and with complex pharmacokinetic proles. These Molecular Targeted Agents (MTA) do not aim at eradicating tumor cells in a short time frame as cytotoxic treatments do, but rather they try to continuously repair the pathological behaviour of tumor cells. Usually, they are administered in a chronic fashion with a high frequency (weekly or even daily). These specicities lead to a change in the nature of related toxic events and therefore, potentially, in the adverse events limiting the dose [START_REF] Booth | On behalf of the Task Force on Methodology for the Development of Innovative Cancer Therapies, Endpoints and other considerations in phase I studies of targeted anticancer therapy: Recommendations from the task force on Methodology for the Development of Innovative Cancer Therapies (MDICT)[END_REF][START_REF] Soria | Phase 1 trials of molecular targeted therapies : Are we evaluating toxicities properly?[END_REF]. In particular, severe late onset toxicities could appear due to drug accumulation following a chronic schedule. Prolonged or repeated moderate toxicity are also relevant for consideration: some moderate toxicities, although deemed to be tolerable when single, signicantly impact patient health when becoming chronic, acquiring an intolerable character [START_REF] Soria | Phase 1 trials of molecular targeted therapies : Are we evaluating toxicities properly?[END_REF][START_REF] Postel-Vinay | Phase I Trials of Molecularly Targeted Agents: Should We Pay More Attention to Late Toxicities?[END_REF]. Therefore, these late onset toxicities and moderate grade repeated adverse events should be considered when MTD is established. When dealing with compounds causing such toxicities, MTD assessment based on binary endpoint of toxicity occurring within a short time period leads to the recommendation of a nal dose regimen which is not tolerated well in the long-term.

Therefore, for MTA, even if early Dose Limiting Toxicities (DLT) are still largely used for driving dose escalation, experts have come to an agreement that, at the end of the study, a dose should be recommended on the basis of any major toxic events whatever the time window within which they appear [START_REF] Postel-Vinay | Phase I Trials of Molecularly Targeted Agents: Should We Pay More Attention to Late Toxicities?[END_REF]. Models are then needed to predict late, repeated or prolonged events with their severity to help in the assessment of the recommended dose, and thus process all information conveyed by the data.

In this paper, we consider the position of data prediction for dose recommendation at the end of study when no more data are awaited. The purpose is not to try to deal with incompleteness of the data as in an ongoing study framework but rather to propose models able to predict any relevant safety data with their time of occurrence as a function of dose for nal dose recommendation purposes.

Dierent approaches have already been proposed for working with late toxicity data, even if time is still needed to translate this research based knowledge into operational practice. Time-to-Event Continuous Reassessment Method (TITE-CRM) [START_REF] Cheung | Sequential designs for phase I clinical trials with late-onset toxicities[END_REF] is probably the simplest of the proposed methods. It was the rst attempt to incorporate information coming from late event in partially observed subjects. This method consists in using a binomial likelihood weighted according to the proportion of observation duration. Another interesting method is the time to event approach. Commonly used to model ecacy data, the time to event approach has also been used for safety data in the context of dose-nding trials [START_REF] Scala | A Novel Bayesian Dose-Escalation Phase Ib Design Investigating Safety of Combination of RAD001 (Everolimus) With Chemotherapy Plus Trastuzumab in Patients With HER2-Overexpressing Metastatic Breast Cancer With Prior Resistance to Trastuzumab[END_REF]. A 2parameter Weibull distribution is used to model the occurrence of DLT at any time during study with a DLT log-hazard function increasing linearly with the dose. Toxicity scores (TS) is another approach which consists in dening and modelling a score describing the intensity of the cumulated toxicity events observed in patients. By setting a target score, MTD is assessed. This approach is derived in dierent ways by considering: 1. dierent toxicity score functions, such as the toxicity index proposed by Rogatko [START_REF] Rogatko | Patient characteristics compete with dose as predictors of acute treatment toxicity in early phase clinical trials[END_REF], the Equivalent Toxicity Level [START_REF] Chen | A novel toxicity scoring system treating toxicity response as a quasi-continuous variable in Phase I clinical trials[END_REF] or the toxicity burden score [START_REF] Lee | Toxicity burden score: a novel approach to summarize multiple toxic eects[END_REF] 2. dierent regression models, such as the probit model [START_REF] Lee | Continual reassessment method with multiple toxicity constraints[END_REF], the Beta regression model [START_REF] Pottho | Flexible Phase I Clinical Trials: Allowing for Nonbinary Toxicity Response and Removal of Other Common Limitations[END_REF], or the Multinomial regression [START_REF] Bekele | Risk-Group-Specic Dose Finding Based on an Average Toxicity Score[END_REF]. TS methods allow an assessment of the toxicity load taking into account the type and the severity of the adverse events regardless of the time they occur, since a toxicity score function can be elicited.

Even if all of these methods are interesting for MTD assessment, only the model proposed by Doussau et. al [START_REF] Doussau | Dose-nding design using mixed-eect proportional odds model for longitudinal graded toxicity data in phase I oncology clinical trials[END_REF] is able to model at the same time, the severity of an adverse event, the time to the occurrence of the event and its duration. TITE-CRM and TS methods do not allow time to event occurrence to be modelled. TITE-CRM approach and time to event methods do not allow the severity of the events to be modelled.

In this work, we propose to explore dierent marginal and Markov models able to predict longitudinal categorical data as is typically observed in early phase trials for the nal assessment of MTD. Both population-averaged and mixed-eect versions of these models are investigated. In the next section, the statistical models under assessment will be described. Section 3 reports the simulation plan allowing the comparison of the competing models and the corresponding results. Finally, we discuss in Section 4 how the retained models can help in the recommendation of a dose at the end of a typical phase I trial.

Statistical models

Let N be the number of patients included in the clinical trial, and let Y ij be the random variable of toxicity observation for individual i = 1, . . . , N at j th time point after drug intake, j = 1, . . . , n i . For the sake of simplicity, the Y ij 's are assumed regularly spaced in time. Let

Y i = (Y i1 , . . . , Y i,n i )
, where n i is the total number of observations for patient i. Note that the number of individual toxicity observations n i is not necessarily the same for all patients, typically when early withdrawals occur during the clinical trial. Most often when dealing with toxicities, Y ij is an ordered categorical variable taking integer values from 0 (no severe toxicity events) to K (most severe toxicity events) corresponding to a grading of the event severity. Implicitly, an initial observation at time 0 (just before drug intake) Y i0 is assumed, but not measured. It will be considered hereinafter that Y i0 ≡ 0 for all i = 1, . . . , N . This implies that immediate toxic reactions at very rst drug administration will only be explained by baseline covariates, since treatment duration and cumulative drug exposure will be both close to zero. Hereafter, such an assumption on Y i0 will be especially useful for Markov approaches for full specication of the corresponding models. Let x ij = (x 1 ij , . . . , x p ij ) ∈ R p denote some possibly time-varying explanatory variables for individual i (such as the time to the beginning of the treatment of the j th record) and let x i = (x ij ) j∈1,...,n i be the matrix of explanatory variables of subject i for the whole observation period with possibly subject-specic parameter ϕ i . The N patients are reasonably assumed to be independent, thus two observation variables Y ij |x i , ϕ i and Y i j |x i , ϕ i would be considered independent when i = i hereinafter. The objective of the present section is to propose statistical models for the Y ij 's which adequately describe and explain the variability in the occurrence of toxic eects. As Figure 1 exemplies, this variability is typically split into some between-patient variability (ie some patients are more susceptible to the toxic eects than others) and some within-patient variability (ie toxicity manifestations for a given patient vary over time). Specically, toxic eects induced by MTA-based treatments turn out to be persistent, suggesting some dependence between the observations of a single patient. The graphs in Figure 1 show that when a toxic eect appears, it is more likely to persist several weeks. The usual statistical models proposed for the analysis of early clinical trials in oncology are not longitudinal data models, and therefore, they do not account for this dependence between the observations of each patient. Doussau et al [START_REF] Doussau | Dose-nding design using mixed-eect proportional odds model for longitudinal graded toxicity data in phase I oncology clinical trials[END_REF] proposed a mixedeect categorical model for toxicity data. Several modelling approaches are suggested in what follows through either population-averaged models or random eects models. For the sake of simplicity, these modelling approaches are presented in detail for the simplest case of binary responses (K = 1) but they could be easily extended to any value of K as briey shown in subsection 2.2.

Models for binary longitudinal data 2.1.1 Population-averaged models

In population-averaged models, the relation between the observations and the subject-specic explanatory variables is assumed to be the same for all subjects. This relation is entirely parametrized by the set of parameters ϕ ∈ R q which is common to the N subjects (this is equivalent to setting ϕ i ≡ ϕ for all individuals i ∈ 1, . . . , N in the population). With such an approach, the variability between subjects is supposed to be fully (and only) described by means of the explanatory variables. Two population averaged modelling strategies are mainly referenced in the literature to model longitudinal data: marginal models and transition models (see for instance [START_REF] Diggle | Analysis of Longitudinal Data[END_REF] and [START_REF] Jansen | Analyzing Incomplete Discrete Longitudinal Clinical Trial Data[END_REF]). The use of either marginal or transition models depends on the main scientic question in hand and corresponds to dierent handlings of the within-subject variability in the data, as we shall see hereafter. Thus, the major dierence between these two modellings lies in the interpretation of the regression coecients.

a) Markov models

Focusing on time transitions, Markov models are a possible way to describe longitudinal data: one is then interested in how some explanatory variables inuence the change of toxicity level over time. The idea of parsimony behind these models is that what will happen in the future is only conditioned by the present state of the system, not by its entire past history. Since some toxic eects might be persistent in time, this seems to be a reasonable assumption in oncology. Here, one is interested in toxic eects induced by a given treatment. Since according to our notations j = 0 matches the beginning of the treatment period, it can reasonably be assumed that treatment-related toxic eects are never recorded at time t i0 , thus the following known initial state for the Markov chain:

P(Y i0 = 0) = 1. (1)
Extensions to any other (possibly ϕ-parametrized) initial distributions would naturally be straightforward. Let us denote H ij = (Y i,j-1 , . . . , Y i1 ) the past recorded trajectory until record Y ij with the convention that H i1 = ∅. Additive time homogeneous transition models could be specied in the following general form:

P(Y ij = 1|H ij , ϕ) = F (µ + p =1 β x ij + f (H ij , α)), (2)
where

f (H ij , α) = min(q,j) r=1 α r 1(Y i,j-r = 0) if j > 1, f (H i1 , α) = 0. (3)
Here, parameters µ, β and α r are unknown constants making up the parameter set ϕ, F stands for any cumulative distribution function, and q is the order of the Markov chain. First-order Markov models (q = 1) are the most widely used transition models. By noting p ϕ ij (1|k ) the probability

P(Y ij = 1|Y i,j-1 = k , ϕ), k = 0, 1, equation (3) becomes p ϕ ij (1|0) = F (µ + p =1 β x ij + α 1 ) , p ϕ ij (1|1) = F (µ + p =1 β x ij ), (4)
the transition matrix from time j -1 to time j of subject i can be easily derived as:

Q ϕ ij = 1 -p ϕ ij (1|0) p ϕ ij (1|0) 1 -p ϕ ij (1|1) p ϕ ij (1|1) (5)
and the full set of (unknown) model parameters is given by ϕ = (µ, β 1 , . . . , β p , α 1 ).

In transition models dened through equation ( 2), the joint distribution of the data vectors Y i can be easily derived through the Markov property that lends a parsimonious form to the general decomposition:

P(Y i |ϕ) = P(Y i0 |ϕ) n i j=1 P(Y ij |H ij , ϕ). (6)
Thus, likelihood-based methods as well as moment-based ones can be implemented to estimate the unknown parameter ϕ from the observations. In particular, the reader's attention is drawn to the fact that here, equation ( 2) species a Generalized Linear Model (GLM) for the conditional distribution of Y ij given the past responses H ij . Indeed, equation ( 2) regresses Y ij on an extended set of explanatory variables which combines the x ij 's and H ij . Thus, it is possible to proceed with estimation of the transition parameters ϕ using inferential techniques specic to GLMs for independent data. Since the resolution of the estimating equations is generally intractable in GLMs, some numerical tricks are required. The well-known Iterative Weighted Least Squares (IWLS) method is such a technique [START_REF] Mccullagh | Generalized Linear Models, Second Edition[END_REF]. It is implemented in the standard glm function of the R software.

b) Marginal models with incomplete specication

When choosing marginal models, by denition, one focuses on the potential eect of explanatory variables on the marginal distributions of the individual temporal observations Y ij (theoretically obtained after integrating out all other Y ij , j = j in equation ( 6)). Marginal models are usually stated in two steps: the marginal distributions of the Y ij 's are modelled separately from the within-subject correlations. In a very general setting, marginal models could be written as follows:

• Step 1: specication of the marginal distribution of each variable Y ij as a function of explanatory variables

P(Y i,j = 1|ϕ) = F (µ + p =1 β x ij ) , 1 ≤ i ≤ N , 1 ≤ j ≤ n i , (7)
•

Step 2: specication of the correlation structure between two toxicity observations Y ij and Y ik from the same subject

Corr(Y i,j , Y i,k |ϕ) = ρ jk (α) , 1 ≤ i ≤ N , 1 ≤ j, k ≤ n i . (8)
Here, F could be any cumulative distribution function, parameters µ and β are constants to be estimated, and ρ jk is a known function of an unknown set of parameters α which species the nature of the autocorrelation. One specic example might be the autoregressive autocorrelation model:

ρ jk (α) = α |j-k| , α ∈ [0, 1], (9)
which assumes that measurements spaced in time are positively correlated, but less correlated than measurements close in time. Many other autocorrelation models can be chosen depending on the user's assumptions concerning the dependence between individual observations. Again, setting β = (β 1 , . . . , β p ), the full set of (unknown) model parameters is given by ϕ = (µ, β, α). Note that, although similar notations are used for the parameters as in the Markovian case, they do not entail the same interpretation (see subsection 2.3 below).

In addition, with the most remarkable exception of independence

P(Y i |ϕ) = n i j=1 P(Y ij |ϕ)
, equations ( 15) and ( 16) are not sucient to specify the entire model properly. Indeed, it is generally not possible to derive the joint distribution of the whole vector of individual observations Y i from ( 15) and ( 16). Thus, likelihood methods cannot be implemented to estimate the unknown parameters ϕ from the observations. For this reason, a well-known alternative comes with the resolution of Generalized Estimating Equations (GEE), see for example [START_REF] Zeger | Longitudinal data analysis for discrete and continuous outcomes[END_REF] for details. Marginal models also fall within the scope of marginal quasi-likelihood and penalized quasi-likelihood methods introduced in [START_REF] Breslow | Approximate Inference in Generalized Linear Mixed Models[END_REF] and [START_REF] Wolnger | Generalized linear mixed models a pseudo-likelihood approach[END_REF] respectively.

Random eects models

In many cases however, the knowledge of the x ij 's is not sucient to fully describe the interpatient variability. In such cases, it is preferable to adopt a subject-specic approach in which the responses (ie the temporal toxicity observations) are modelled as a function of covariates and subject-specic random eects ϕ i . Such models are usually dened hierarchically by dening the ϕ i 's as random variables, so that population-averaged conclusions can still be derived through the so-called population parameters θ. Dening the individual parameters as random variables is also a way to handle the correlations between temporal observations from each subject. Random-eects models are specied in two stages:

• Stage 1 : the conditional distribution of the individual observations Y i given the individual parameters ϕ i assumes the same kind of relation between the temporal observations and the explanatory variables for the N subjects:

Y i |ϕ i , x i ∼ p(•|ϕ i , x i ), (10)
• Stage 2 : the individual parameters ϕ i are dened as independent random variables, the distribution of which is parametrized by a set of population parameters θ:

ϕ i ∼ i.i.d. p(•|θ). (11)
In most cases, the individual parameters are assumed to be Gaussian random variables. Multiple variations of the random-eects model could be derived for longitudinal categorical toxicity observations. Two specic examples of interest among so many others are given here:

1. The generalized linear mixed-eects model:

P(Y i,j = 1|ϕ i ) = F (µ i + p =1 β i, x ij ) , P(Y i,j , Y i,j |ϕ i ) = P(Y i,j |ϕ i ) × P(Y i,j |ϕ i ) , ϕ i = (µ i , β i,1 , . . . , β i,p ) ∼ i.i.d. N (ϕ, Ω). ( 12
)
2. The mixed-eects rst-order Markov model, an extension of the above transition models

(2) -( 3) in a random-eects approach:

P(Y ij = 1|Y i,j-1 , ϕ i ) = F (µ i + p =1 β i, x ij + α i 1(Y i,j-1 = 0)), ϕ i = (µ i , α i , β i1 , . . . , β ip ) ∼ i.i.d. N (ϕ, Ω). (13)
In both examples, the population parameters are taken as the mean and the covariance matrix of normal population distribution θ = (ϕ, Ω).

Inference of the population parameters θ in mixed-eects models is classically performed by maximization of the model likelihood. Owing to their hierarchical denition through equations ( 10) and [START_REF] Scala | A Novel Bayesian Dose-Escalation Phase Ib Design Investigating Safety of Combination of RAD001 (Everolimus) With Chemotherapy Plus Trastuzumab in Patients With HER2-Overexpressing Metastatic Breast Cancer With Prior Resistance to Trastuzumab[END_REF], the joint marginal likelihood of each individual sequence of observations is straightforwardly given by integrating out p(•, x i |ϕ i ) over the distribution of the random eects:

P(Y i |θ) = p(Y i |ϕ i , x i )p(ϕ i |θ)dϕ i . (14)
In many cases, equation ( 14) does not have any explicit expression. Thus, maximum likelihood is often performed based on linearizations of ( 14) (see for instance chapter 6 in [START_REF] Davidian | Nonlinear Models for Repeated Measurement Data[END_REF]), numerical approximations of the integral or via EM-based algorithms [START_REF] Xu | Likelihood Analysis of Multivariate Probit Models Using a Parameter Expanded MCEM Algorithm[END_REF]. Marginal quasi-likelihood [START_REF] Breslow | Approximate Inference in Generalized Linear Mixed Models[END_REF] and penalized pseudo-likelihood methods [START_REF] Wolnger | Generalized linear mixed models a pseudo-likelihood approach[END_REF] have also been specically developed to perform parameter estimation in mixed-eects generalized linear models. As in the xed-eects setting, by denition of the Markov models given above, mixed-eects Markov models fall within the scope of this method. To t generalized linear mixed models, one could use the function glmer of the R software, which approximates the integral in ( 14) by quadrature methods.

Models for ordinal longitudinal data

In practice, toxicity data may fall into more categories than just two possibilities, the possible toxicity values establishing a ranking between symptoms (see Figure 1 for instance). Some Markov and marginal models for longitudinal ordinal data could be easily derived from the longitudinal models for binary data presented in Section 2.1 whether in a xed-eects setting or in a random eects setting. Without loss of generality, it is considered here that the grades of toxicity are integer values from 0 (no severe toxic event) to K (the most severe toxic events). Some possible models for this setting are listed below, but not detailed. Indeed, their interpretation and inference are basically similar to those for the corresponding binary data models.

(a) Fixed-eects (population averaged) Markov models:

P(Y ij ≤ k|H ij , ϕ) = F (µ k + p =1 β x ij + f (H ij , α)) 0 ≤ k ≤ K -1, P(Y ij ≤ K|H ij , ϕ) = 1, f (H ij , α) = min(q,j) r=1 K-1 k=0 α rk 1(Y i,j-r = k) if j > 1, f (H i1 , α) = 0 otherwise,
where parameters µ 0 < µ 1 < . . . < µ K-1 , β and α rk are unknown constants and compose the parameter set ϕ, F stands for any cumulative distribution function, and q is the order of the Markov chain. As above, for q = 1 (rst-order Markov models), the transition matrix from time j -1 to time j of subject i can be easily derived from the probabilities

P(Y ij ≤ k|H ij , ϕ).
(b) Fixed-eects (population-averaged) marginal models with incomplete specication: As above, marginal models for longitudinal ordinal data are dened in two steps:

•

Step 1: specication of the marginal distribution of each variable Y ij as a function of explanatory variables

P(Y i,j ≤ k|ϕ) = F (µ k + p =1 β x ij ) , 0 ≤ k ≤ K -1, P(Y i,j ≤ K|ϕ) = 1 , 1 ≤ i ≤ N , 1 ≤ j ≤ n i , (15)
• Step 2: specication of the correlation structure between two toxicity observations Y ij and Y ij from the same subject

Corr(Y i,j , Y i,j |ϕ) = ρ jj (α) , 1 ≤ i ≤ N , 1 ≤ j, j ≤ n i . (16)
F stands for any cumulative distribution function, ϕ = (µ 0 , . . . , µ K-1 , β 1 , . . . , β p , α) is the unknown vector of model parameters to be estimated, such that µ 0 < µ 1 < . . . < µ K-1 , and ρ jj is a known function of an unknown set of parameters α which species the nature of the autocorrelation, for example an autoregressive one.

(c) Generalized linear mixed-eects models:

P(Y i,j ≤ k|ϕ i ) = F (µ ik + p =1 β i, x ij ) , 0 ≤ k ≤ K -1, P(Y i,j ≤ K|ϕ i ) = 1 , P(Y i,j , Y i,j |ϕ i ) = P(Y i,j |ϕ i ) × P(Y i, j |ϕ i ) , ϕ i = (µ i0 , . . . , µ i,K-1 , β i1 , . . . , β ip ) ∼ i.i.d. N (ϕ, Ω), ( 17
)
where F is a given cumulative distribution function and θ = (ϕ, Ω) is the set of population parameters.

(d) Mixed-eects Markov models:

P(Y ij ≤ k|H ij , ϕ i ) = F (µ ik + p =1 β i, x ij + f (H ij , α i )) , 0 ≤ k ≤ K -1, P(Y ij ≤ K|H ij , ϕ i ) = 1, f (H ij , α i ) = j-1 r=1 K-1 k=0 α irk 1(Y i,j-r = k), ϕ i = (µ i0 , . . . , µ i,K-1 , β i1 , . . . , β ip , α i ) ∼ i.i.d. N (ϕ, Ω), (18)
where, as above, F is a given cumulative distribution function and θ = (ϕ, Ω) is the set of population parameters.

Some examples

To cast light on the dierence between the previous model families (Markov and marginal models whether in a population-averaged or a random-eects setting), as well as on the nature of this dierence, a few model examples are presented and compared. For the sake of simplicity, the toxicity observations are considered as binary outcomes. This would correspond to cases where the user's interest would be for instance the presence (Y ij = 1) or absence (Y ij = 0) of toxic eects over time or correspond to the appearance with time of disabling (Y ij = 1)/nondisabling (Y ij = 0) toxic eects. Once again, in practice, there may be more categories than just two possibilities, but the binary situation will capture the essence of the example, whilst avoiding arithmetic and notational complexities. For the purpose of the example here, only one time non-varying explanatory variable c i is considered. The main focus will be on the meaning of the regression coecient β related to covariate c i . For that purpose, the expression of the marginal distribution of the temporal observations with respect to the model parameters in each example -entirely given by the probability P(Y ij = 1) since the data are binary -will be introduced. The probit function, denoted as Φ, is chosen for specifying the model examples for the comparative study. Consider

• a) the generalized linear model (regardless of the correlation structure between the observations of each subject) given by

P(Y ij = 1) = Φ (µ + βc i ) ,
• b) the mixed eects generalized linear model such that

P(Y ij = 1|η i ) = Φ (µ + (β + η i )c i ) , η i ∼ i.i.d.
N (0, ω 2 ),

• and c) the Markov model given by

P(Y ij = 1|Y i,j-1 = k) = Φ (µ + βc i + αk) .
In the last two examples, the marginal distribution of Y ij can be derived, although not necessarily explicitly. Moreover, it generally does not lead to an equality with Φ (µ + βc i ) as in the simple GLM case. In the mixed-eects GLM, it requires integration over the random eect η i distribution:

P(Y ij = 1) = Φ (µ + (β + η i )c i ) p(η i |ω 2 )dη i , = Φ (µ + βc i ) ,
whereas in the Markov model, it requires summing over all possible values for the past observations from time 0 up to time j -1:

P(Y ij = 1) = 1 k 0 =0 1 k 1 =0
. . .

1 k j-1 =0 P(Y i0 = k 0 ) j l=1 P(Y il = k l |Y i,l-1 = k l-1 ), = Φ (µ + βc i ) , with P(Y il = k l |Y i,l-1 = k l-1 )
given by Φ (µ + βc i + αk l-1 ) if k l = 1 and 1 -Φ (µ + βc i + αk l-1 ) otherwise.

Figure 2 illustrates how P(Y ij = 1) changes according to the value of β in cases a), b) and c). The larger ω, the more distant the marginal distributions of b) and a) will be. Moreover, while the marginal distribution remains common to all the Y ij 's of a given subject with xed and random eects GLMs, it does dier from time point to time point in the Markov model. This corroborates the fact that the interpretation of parameter β is completely dierent according to the model. In the xed-eects GLM, β characterizes the eect of covariate c i on the probability for observing presence Y i,j = 1 at time j at the population level, whereas in the random-eects GLM, it has a similar meaning at the individual level, i.e. it is conditional to the individual parameters. Population and individual conclusions are of course combined when the betweensubject variability is null. On the contrary, in model c), β is to be interpreted in terms of transitions.

3 Simulation study

Objectives and simulation settings

The simulation study is performed to compare several longitudinal models with a specic focus on their predictive characteristics, on a whole cohort of patients (rst simulation study) as well as on small groups of patients treated with similar dose levels (second simulation study). In both cases, the robustness to model misspecication, thanks to dierent simulation scenarios specied thereafter, is also explored. The whole study is based on original datasets of two similar Phase I studies on an MTA where the toxicity is measured on scale from 0 to 5. Since the models, data generation procedures and predictive criteria are common to the dierent steps of the simulation study, they are described once for clarity before summarizing the results of the two simulation studies.

Model construction from real clinical data

The subsequent simulation studies are based on four models which correspond each to one modelling approach introduced in Section 2. Each one assumes the subjects to be independent of one another. Two of these four models are Generalized Linear Models, with independent and correlated data respectively, and the others two are transitional models, with xed and random eects respectively. These four models are built using two real datasets from similar early clinical trials investigating the same MTA. In these trials, the toxicity observations are binary outcomes (K = 1): Y ij = 0 if the observed toxicity level is 0 or 1, Y ij = 1 otherwise, and the observations are collected at regular times, ie t ij = j for all i = 1, . . . , N and all j = 1, . . . , n i . Many covariates are available to predict the toxicity responses. Some of these covariates are highly correlated as, for example, the dose level (ie the dose administered at each administration) and the dose intensity (i.e. the average administered dose according to time of treatment). Since highly correlated covariates may provide poor parameter estimations and lead to possible misinterpretations, the available covariates are grouped into 4 categories and nal covariate selection for each model is performed by Information criteria (AIC) under restrictions that no more than one covariate from each group may be included in the model. The following covariate categories are considered: -1. Individual characteristics like age or gender -2. Treatment period and time from the last drug administration (TFLDA) -3. Measure of administered treatment like dose level (DL) or dose intensity (DI) -4. Cumulated treatment exposure like cumulated dose (CumDose) or cumulated AUC of patients' pharmacokinetic curves (CumExp). For most of these covariates, the logarithmic transformations have also been proposed (and assigned to their corresponding categories). The following models are selected this way, with Φ denoting for the probit link:

• GLM (Generalized linear model): for each subject i, the temporal toxicity observations are given by:

P(Y ij = 1|ϕ 0 ) = Φ(µ (0) + β (0) 1 log t ij + β (0) 2 TFLDA ij + β (0) 3 log DI i + β (0) 4 CumDose ij ), j = 1, . . . , n i , where ϕ 0 = (µ (0) , β (0) 1 , β (0) 2 , β (0) 3 , β (0) 4 ).
1. First, the ability for each model to predict which patients may or may not have a toxicity event over their period of participation in the trial will be examined. This will be evaluated by means of the Brier score [START_REF] Brier | Verication of forecasts expressed in terms of probability[END_REF], which is a proper score generally used to quantify the quality of such binary predictions. This score reaches its minimum when a perfect matching between the probabilistic forecast and the distribution of the event to be predicted is obtained.

2. The second examined criterion is how accurately each model predict the number of adverse events occurring for each patient during the trial. Since the number of reported toxicities induced by MTAs is likely to be larger for long periods than for short ones, the number of adverse events is scaled according to the duration of the follow-up. The predictions of the number of adverse events are then evaluated with the Continuous Rank Probability Score (CRPS) [START_REF] Epstein | A Scoring System for Probability Forecasts of Ranked Categories[END_REF], which is also a proper score. Thus, the better the model predicts the number of adverse events, the smaller the CRPS will be.

3. Third, the predicted number of weeks with a toxicity response is controlled. As for previous criteria, the number of weeks with a toxicity response is brought back to the duration of the follow-up and the CRPS is used to evaluate the ability for each model to provide such good predictions.

Simulation study 1: comparison on overall predictive characteristics

This rst simulation study aims at comparing the predictive properties of models GLM, AR, F.markov and M.markov without any distinction between the patient dose levels of the MTA of interest. to the prediction of the duration of adverse events and to the prediction of the total number of adverse events respectively. On each gure, the rst line displays the model performance under the three simulation scenarios in cases of small datasets (N = 30) and the second line displays the model performance under the three simulation scenarios in cases of larger datasets (N = 50).

Data generation

The rst noticeable result is that under the bootstrap simulation scenario (Scenario 1)

and whatever the sample size, models AR and M.markov lead to Brier scores and CRPSs smaller than those given by models GLM and F.markov, without much dierence between the performance of AR and M.markov. This especially means that the assumption on the nature of the dependence between longitudinal toxicity observations underlying these two models provides an adequate description of the toxic reactions.

Simulation scenarios Scenario 2 and Scenario 3 are aimed at evaluating the robustness of the candidate models to model misspecication. Under Scenario 2, data are generated according a mixed-eects Markov model. Figures 3, 4 and5 show that the best predictive model according to the three chosen criteria is either the xed-eects Markov model (F.markov) or the mixedeects Markov model (M.markov) when N = 30 and always the mixed-eects Markov model when N = 50. The diculty of recovering the true model from small datasets is probably linked to the diculty of correctly evaluating the between-patient variability in this case but the main assumption of a Markov structure of dependence between observations is sucient to provide correct predictions. Under Scenario 3, data are generated according to a marginal model with an autoregressive correlation structure. Under this scenario, no clear dierence is evidenced between the four models whatever the sample size. Indeed, the median and the dispersion of the three scores are quite similar from one model to another.

As an overall conclusion, this rst simulation study shows that assuming a dependence structure (autoregressive or Markov) between longitudinal toxicity records leads to the most appropriate predictions of the toxic reactions to MTAs. More precisely, Markov models, in particular mixed eects Markov models, would be a robust model choice for accurately describing longitudinal toxicity records. corresponding theoretical probabilities of having at least one toxic event over the 15 weeks of monitoring are 0.06, 0.11, 0.19, 0.3, 0.41 respectively for doses DI 1 , DI 2 , DI 3 , DI 4 and DI 5 under Scenario 4 and 0.08, 0.19, 0.32, 0.42, 0.50 under Scenario 5.

Each pair of simulated datasets includes, as above, one training set and one validation set.

Whatever the scenario used for data simulation (Scenario 4 or Scenario 5), the procedure is the following. For each of the 1000 simulated training cohorts, the two competing models (AR and M.markov) are rst inferred; then, using the estimated model parameters, the predicted probability of having at least one toxic event is computed for each patient of the corresponding validation cohort. The empirical distributions of the predicted probabilities of having at least one toxic event over the 15 weeks of monitoring can be deduced from the 1000 simulated cohorts. It is then compared to the theoretical probabilities stated above using the CRPS.

Results

The results are depicted Figure 6. The rst line gives the results obtained when simulating balanced dose level groups, whereas the second line displays the results arising from simulating unbalanced dose level groups. On the left, data are simulated according to a mixed-eects hidden Markov model (Scenario 4 ); on the right, data are simulated according to a generalized linear model with an autoregressive correlation structure (Scenario 5 ). In case a) (balanced dose level groups), model AR provides better predictions on low dose intensities (DI 1 and DI 2 ) whatever the scenario used for data simulation, while model M.markov tends to provide better predictions on higher dose intensities (from DI 3 to DI 5 ), which are closer to the MTD (e.g. 35%). In case b) (unbalanced dose level groups), model M.markov provides better overall predictions whatever the scenario used for data simulation. Indeed, the CRPS distributions for each dose group show a lower median score and smaller variability in prediction when obtained by inferring the autoregressive model.

Discussion

Categorical and temporal extensions

One can imagine many extensions of the longitudinal data models described above. First of all, although the illustrative examples presented in Sections 2.3 and 3.1.1 restrict for simplicity the temporal toxicity responses to binary outcomes (typically presence or absence of side eects), we emphasize that the approach also applies to more general types of outcomes like ordinal categorical ones. The model can consider dierent grades of toxicity, for instance from 0 to 5 according to NCI-CTCAE nomenclature, which would allow for a more sophisticated description of the evolution of toxicity over time with respect to the clinician's interest. Secondly, our presentation focused on modelling some unidimensional temporal toxicity responses. However, just like standard oncology therapies, MTAs are likely to induce several kinds of toxic eects like, without being exhaustive, ocular toxicities, peripheral neuropathies or urinary toxicities. When determining the Maximum Tolerated Dose, one should take into account the way all these toxicities jointly occur over time, for which the unidimensional longitudinal models presented in Section 2 oer limited prospects. There is indeed no reason to consider that the underlying mechanisms for these dierent categories of toxicities are the same or that these toxic eects arise independently. Therefore, it is of strong interest to ultimately extend the models above to (possibly correlated) multidimensional temporal responses. Data augmentation as in [START_REF] Albert | Bayesian Analysis of Binary and Polychotomous Response Data[END_REF] is the convenient entry point to develop hierarchical stochastic structures with increased complexity: introducing a latent multivariate (normal autoregressive) layer is the price to pay in order to match all these operational requirements.

By nature, the serial dependence is implemented in the Markov models we considered, in addition to a between-individual variability induced by the random eects for the last one. A hidden Markov model (HMM) may be an attractive alternative for the case study ( [START_REF] Cappe | Inference in Hidden Markov Models[END_REF]). Instead of establishing a direct transition between observable states Y ik as given by equation 2, in HMM, the time dependence stems from a latent Markov model Z. For instance a rst order Normal AR model would be written:

Z ij = ρ i Z i (j-1) + σ i ε ij ε ij ∼ N (0, 1).
Although ρ i can be high to keep the latent process Z ij tight enough during long sequences, some dispersion is nevertheless allowed thanks to the stochastic observation equation given by:

P(Y ij = 1|Z ij ) = Φ (µ + βc i + Z ij )
Such HMM with random eects ((ρ i , σ i ) ∼ p(.|θ)) have recently been proposed by [START_REF] Altman | Mixed hidden Markov models: An extension of the hidden Markov model to the longitudinal data setting[END_REF] in the advanced statistical literature.

On the path of personalized medicine

The present study mainly focused on the interests of longitudinal data models at the population scale, the idea being to depict the mechanism of the appearance of toxic eects over time and their variations around some population average. But the other aspect of interest for longitudinal data models lies at the individual level. Once the model is inferred from the whole trial data, one can take advantage of its predictions at the individual levels. Borrowing strength from neighbours but allowing for some heterogeneity among individuals, random eect models seem here to be ecient tools purposely tailored.

Missingness

Getting personalized knowledge about individuals appears particularly relevant in the context of phase I clinical trials in oncology since overly strong adverse events lead to early withdrawal from the trial. It is worth noting that such early drop-out cannot be considered independently of the observable and missing data. Describing the missingness mechanism yields additional modelling and inference challenges not dealt with this paper as well as extra programming not available in standard software. [START_REF] Jansen | Analyzing Incomplete Discrete Longitudinal Clinical Trial Data[END_REF] provides an overview of the various modelling frameworks for non-Gaussian longitudinal data, with non ignorable missingness processes. Once a patient shows a critical toxicity prole, the model can help predict his near future and anticipate stronger adverse eects and can therefore possibly adapting his treatment for more comfort.

Words of caution

The introduction of longitudinal models in contexts of phase I clinical trials in oncology however is not exempt from some caveats. First, longitudinal models allow the use of time-varying covariates to predict the toxicity outcomes, which was not possible with approaches based on a one-shot outcome (non longitudinal models). Clearly, these time-varying covariates are a very rich source of information but also a potential source of confusion with the temporal memory of toxicity responses. It is however up to the user to use these time-varying covariates cautiously in a way to avoid misinterpretation of the model parameter estimations. A typical situation where some problems may arise is when the treatment is readjusted (delayed administration or reduced dose administration) because of the occurrence of disabling or persistent toxic eects. In such a situation, the causal relationship between toxicities and drug exposure becomes ambiguous since it works in both directions.

A second warning concerns the amount of data required to satisfactorily infer complex longitudinal models. The more sophisticated the model becomes, the more data is needed to get trustworthy parameter estimations. Yet the size of studied groups in phase I clinical trials is often rather limited, from say 20 to 80 patients. As usual, model complexity has to be weighted with sample size.

Conclusion

From retrospective analysis to experimental design

In this paper, the retrospective analyses based on data collected from phase I clinical trials in oncology plead in favour of the use of longitudinal models, especially Markovian models with random eects. We believe that these advanced stochastic structures are not limited to the statistical description of the data, but in this section we show that they could also develop as operational tools. Our work was particularly motivated by the recent development of MTAs, which rely on more complex mechanisms of action than do standard oncology treatments. Multiple dose administrations to the same patient are performed over multiple cycles. MTAs provoke repeated side eects of variable intensity in the long term that are recorded in the follow-up of each patient (for instance such follow-up may last over six cycles of one month). In the practice of clinical trials however, the next cohort of patients to be treated is launched long before all the data of the previous cohort has been completely recorded, say on a monthly basis. To cope with this short horizon, most studies currently reduce the data from the previous cohorts to a binary end point, the occurrence of a toxicity, ignoring the incompleteness of the data, the recurrence pattern of events and the various grades of toxicity. Such preliminary "syntheses" to sum up the follow-up data yield a substantial loss of information, particularly on the temporal evolution of toxicities, and therefore present the risk of recommending an unadapted dose. Longitudinal data models give clinicians the opportunity to take advantage of all available information from the clinical trial. Following the same avenue of thought of the continuous reassessment method, more rened dose-nding strategies could be elaborated after updating the knowledge about the longitudinal model parameters as information accumulates and late toxicities from previous cohorts come out.

Avoiding dose limiting toxicities

Up to now, unless a dose limiting toxicity occurs, the dose received by a patient is usually not changed from one cycle to the next. Longitudinal models intend to assess the mechanisms underlying the appearance of toxicities over time and allow for the incorporation of timevarying covariates. They account for the dependence between consecutive observations and thus make us able to assess dierent toxicity proles like late, repeated or persistent toxicities. We show that these models also have interesting predictive properties for separate pre-dened dose categories. Therefore the extra information gained from each cycle provides more precise estimation to understand the individual patient's behaviour with respect to toxicity during the follow-up period and might even enable a better adaptation of the dose within the cycles of a given patient.

Redening the MTD

However over-simplistic the statistical approaches presently retained for phase I clinical trials in oncology may be, they have the great merit of designing clinical trials targeting MTD. Nevertheless, as a direct end-product of these simple models, the MTD is in essence a statistical quantile. The single summarizing outcome per patient follow-up is considered as a Bernoulli event, and the MTD is merely the dose corresponding to the 33% quantile of the probability model that rules such Bernoulli events, with one shot per patient. Now, phase I oncology experts have come to the agreement that new recommendations for toxicity evaluation are needed, in particular requiring the assessment of new MTD based on multiple endpoints [START_REF] Soria | Phase 1 trials of molecular targeted therapies : Are we evaluating toxicities properly?[END_REF][START_REF] Postel-Vinay | Towards new methods for the determination of dose limiting toxicities and the assessment of the recommended dose for further studies of molecularly targeted agents -Dose-Limiting Toxicity and Toxicity Assessment Recommendation Group for Early Trials of Targeted therapies, an European Organisation for Research and Treatment of Cancer-led study[END_REF].

In the present work, we study alternative models allowing us to learn from multivariate longitudinal categorical data. Unlike previous one-shot Bernoulli models, such longitudinal models oer the appropriate framework to enrich the denition of the Maximum Tolerated Dose (which remains the main objective of a Phase I clinical trial in oncology) by taking into consideration dierent endpoints longitudinally followed. A renewed denition of the Maximum Tolerated Dose appropriate to MTA compounds should take into consideration acute, late-onset, repeated or prolonged toxicities of many kinds along with their related severity. For instance, the MTD could be seen as the maximum dose verifying that dierent toxic events occur with some predened controlled risks:

arg d P(E 1 |d) < t 1 ∧ P(E 2 |d) < t 2 ∧, ..., ∧ P(E S |d) < t S
where E s (s = 1 . . . S) are the toxic events of interest and t s the consented risks of occurrence of these events. As an example, E s could be a late onset peripheral neuropathy event with a grade above or equal to 3 according to NCI-CTCAE, or a Grade 2 to 4 ocular toxicity event lasting more than 4 weeks, or even the simultaneity of those two toxic events. Inasmuch as enough data becomes available to perform a relevant inference, the proposed longitudinal data models allow MTD to be assessed based on the richer denitions that phase I oncologists are willing to promote.

Figures 3 ,

 3 Figures 3,[START_REF] Ji | A modied toxicity probability interval method for dosending trials[END_REF] and 5 display the results related to the prediction of at least one toxic event, to the prediction of the duration of adverse events and to the prediction of the total number of adverse events respectively. On each gure, the rst line displays the model performance under the three simulation scenarios in cases of small datasets (N = 30) and the second line displays the model performance under the three simulation scenarios in cases of larger datasets (N = 50).

Figure 1 :

 1 Figure1: Typical toxicity records from an early phase clinical trial of a specic MTA. Each graph corresponds to a patient and represents the evolution over weeks of the intensity of toxic eects after drug intake on a scale from 0 (no toxic eect) to 4 (disabling eect).

Figure 2 :Figure 4 :Figure 5 :

 245 Figure 2: Relation between P(Y ij = 1) and β in models a) (red), b) (blue) for dierent values of ω (the greater ω is, the more the thicker the line is) and c) (black dotted lines) for dierent values of α. Each graph corresponds to a dierent time point j (j = 1 on the left, j = 2 on the right). The following values are used: µ = 0.1 and c i = 1.

  Three simulation scenarios are considered. To simulate realistic data, individual proles (dose level, drug exposure, time of treatment, etc) are sampled with replacement from the original datasets. Under Scenario1, a sort of bootstrap scenario, the patients' proles (i.e. values for the covariates and length of follow-up) and the corresponding toxicity longitudinal responses are sampled with replacement from the original dataset. Such a simulation technique will help identify the best predictive model without any idea of the true model. Two other simulation scenarios, Scenario2 and Scenario3, are then used to investigate the robustness to model misspecication. In both, the simulated patients' clinical characteristics are obtained by random sampling with replacement from the initial data as in Scenario1, but the corresponding toxic-

ity longitudinal responses are simulated according to mixed-eects Markov model M.markov

Acknowledgements The authors thank Joan Sobota for her advice and comments.

• AR (Marginal model with autoregressive correlation structure): for each subject i, the moments of the temporal toxicity observations are assumed to be such that: P(Y ij = 1|ϕ 1 ) = Φ(µ (1) + β

(1)

4 log CumExp ij ), Corr(Y i,j , Y i,k |ϕ 1 ) = ρ |k-j| , 1 ≤ j, k ≤ n i , where ϕ 1 = (µ (1) , β

2 , β

4 , ρ).

• F.markov (Fixed-eects Markov model): for each subject i, the temporal toxicity observations are given by:

2 log DI i + β

(2)

3 CumExp ij + α (2) Y i,j-1 ),

where ϕ 2 = (µ (2) , β

1 , β

3 , α (2) ).

• M.markov (Mixed-eects Markov model): given some subject-specic parameters ϕ i , the temporal observations of each subject are given by the same structural model:

i + β

(3)

i1 Y i,j-1 ).

Here, the individual parameters are ϕ i = (µ

i1 , β

i2 , β

i3 , β

i4 , α

i1 ). To model the between-patient variability, the simplifying assumption that only parameter α

α ), whilst the other individual parameters are shared among subjects with ignorable variation such that ∀iµ

≡ β

(3) , = 1, . . . , 4.

Thus, the population parameters are θ = (µ (3) , (β

) =1,...,4 , α

1 , ω 2 α ).

Parameter estimation in GLM and F.markov is provided by maximization of the likelihood (Iteratively Re-weighted Least Squares algorithm [START_REF] Mccullagh | Generalized Linear Models, Second Edition[END_REF]), AR parameters estimation is based on a Penalized Quasi-Likelihood maximization method [START_REF] Wolnger | Generalized linear mixed models a pseudo-likelihood approach[END_REF][START_REF] Breslow | Approximate Inference in Generalized Linear Mixed Models[END_REF] and in M.markov, it is given by maximization of the likelihood through the Gauss-Hermite quadrature method [START_REF] Pinheiro | Approximations to the Log-Likelihood Function in the Nonlinear Mixed-Eects Model[END_REF].

Predictive criteria

The ensemble (i.e. whatever the dose level) predictive performance and the dose-specic predictive performance of models GLM, AR, F.markov and M.markov will be assessed according to the three criteria below.

(using the parameters values µ (3) = -6.85, β

and ω α = 0.48) and marginal model with autoregressive correlation structure AR respectively (using the parameters values µ (1) =, β

4 = and ρ =). Since marginal models are rarely completely specied, it is hard to simulate data according to such models. Thus, the simulation technique introduced in [START_REF] Biswas | Generating correlated ordinal categorical random samples[END_REF] is used to generate the data under Scenario3.

For each of these three data generation scenarios, the procedure is the following. 1000 pairs of independent datasets with N patients each are simulated. The rst one (T r k ) is a "training sample" and the second one (T e k ) is a "test sample", k = 1, . . . , 1000. For each value of k, the four models GLM, AR, F.markov and M.markov are inferred from T r k regardless of the scenario used for generating the data, whereas the predictive properties of the models are computed by Monte Carlo from T e k by comparing each patient's responses from T e k with 1000 probable longitudinal toxicity responses under the estimated model. Two sample sizes are investigated: N = 30, which corresponds more or less to the size of most phase I studies, and N = 50, which is the size of the original datasets used to design the present study.

Results

Simulation study 2: comparison on dose-specic predictions

Beyond the comparison of models AR and M.markov according to their overall performances in predicting toxic events for entire cohorts of patients, their predictive performance is also compared on the scale of subgroups of patients having been given the same dose levels. This second simulation study aims at assessing the ability of each model to discriminate between doses, thus at assessing the appropriateness of each model in dose-nding.

Data generation

Two situations are investigated: a) an ideal case where the dose groups are balanced, and b) a more realistic case where the patients' allocation to dose groups is not balanced. 5 dose groups are considered, further denoted as DI 1 , DI 2 , DI 3 , DI 4 and DI 5 . In situation b), the number of simulated patients in each dose group varies from one simulated dataset to another, since in this case, dose level is assigned by randomly choosing a patient and his dose level among those of the original datasets. In both balanced and unbalanced congurations, 1000 pairs of cohorts of 50 patients are simulated assuming the same monitoring duration for all patients, including 9 weeks of treatment and 6 weeks of follow-up. The patients' drug exposure (i.e. the values over time of explanatory variable CumExp) is simulated according to a loglinear model, previously inferred from the original datasets:

96 and ε ij ∼ N (0, 0.16 2 ). Their dose intensities are also simulated. As in previous simulation study, robustness of predictions by dose level to model misspecication is assessed by simulating the patients' toxicity events according to dierent scenarios. Under Scenario 4, a mixed eects Markov model (M.markov) is used to simulate the toxicities, with the following parameter values:

= 0.025, ω 2 α = 1; and under Scenario 5 the toxicities are simulated according to a marginal model with autoregressive correlation structure (AR), using parameter values µ (1) = -1.78, β

= -0.128, β