- Email: [email protected]
A tampered Brownian motion process model for partial step-stress accelerated life testing
Journal of Statistical Planning and Inference 94 (2001) 15–24
www.elsevier.com/locate/jspi
A tampered Brownian motion process model for partial step-stress accelerated life testing a Division
Yili Lua; ∗ , Barry Storerb of Statistical and Mathematical Sciences, Lilly Research Laboratories, Eli Lilly and Company, Lilly Corporate Center, Indianapolis, IN 46285, USA b Fred Hutchinson Cancer Research Center, USA
Received 28 July 1998; received in revised form 15 August 2000; accepted 29 August 2000
Abstract In partial step-stress accelerated life testing, models extrapolating data obtained under more severe conditions to infer the lifetime distribution under normal use conditions are needed. Bhattacharyya (Invited paper for 46th session of the ISI, 1987) proposed a tampered Brownian motion process model and later derived the probability distribution from a decay process perspective without linear assumption. In this paper, the model is described and the features of the failure time distribution are discussed. The maximum likelihood estimates of the parameters in the model and their asymptotic properties are presented. An application of models for step-stress accelerated life test to elds other than engineering is described and illustrated by applying the c 2001 Elsevier tampered Brownian motion process model to data taken from a clinical trial. Science B.V. All rights reserved. MSC: 60J65; 62N05 Keywords: Partial acceleration; Brownian motion process; Inverse Gaussian distribution; Maximum likelihood estimation
1. Introduction Accelerated life testing of material is used to collect failure data quickly when the lifetime of a specimen under normal use conditions is too long. The assumption for using an overstress accelerated life test (ALT) is that the lifetime is adversely a ected by some physical conditions which are referred to as stresses. Accelerated life tests which involve higher levels of physical stress than would be experienced under normal conditions will hasten the process of failure. Step-stress ALT is a special kind of ALT in which the stress on each specimen is increased step-by-step over time. When a test ∗
Corresponding author. Tel.: +1 317-277-3954; fax: +1 317-276-5064 E-mail address: lu [email protected] (Y. Lu).
c 2001 Elsevier Science B.V. All rights reserved. 0378-3758/01/$ - see front matter PII: S 0 3 7 8 - 3 7 5 8 ( 0 0 ) 0 0 2 1 5 - 9
16
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
involves two levels of stress with the rst one as the normal one and has a xed time point for changing stress, it is referred to as a partial step-stress ALT. This will be the setting that we deal with throughout this paper, hence, the word partial is omitted hereafter. The goal of statistical inference for step-stress ALT is to estimate the failure behavior of the specimen in conditions of normal use using the failure data obtained under more severe conditions. To achieve this goal, proper statistical models are needed. The pioneering work for step-stress ALT modeling is by DeGroot and Goel (1979) and Nelson (1980). The former proposed a tampered random variable model and the latter a cumulative exposure (CE) model. Alternatively, Bhattacharyya and Soejoeti (1989) suggested a tampered failure rate model, and Doksum and Hoyland (1991) brought forth a model from the decay process point of view. In Doksum and Hoyland’s approach, instead of modeling the failure time distribution directly, a damage growth process is modeled in relation to the acceleration of the stress. The failure time distribution is then obtained as the distribution of the rst passage-time of the process with respect to a critical boundary. In the meantime, nonparametric models have also been proposed (Shaked and Singpurwalla, 1983; Schmoyer, 1991). Most of the parametric models listed above share a common mathematical assumption that there is a linear association between the failure time before and after a change of the stress. For instance, it is assumed that the decay process under the higher stress is linearly associated with the one under the lower stress in Doksum and Hoyland’s model, which yields a tampered random variable model for the failure time with an inverse Gaussian distribution. It can also be shown that the tampered random variable models compose a proper subset of the cumulative exposure (CE) models, when the CE models are limited to those within location-scale families. Therefore, the availability of non-linear models is very limited when the linear assumption fails in the real world. This paper examines another approach for modeling step-stress ALT. The idea of the approach was mentioned by Bhattacharyya (1987) and the model was later derived by the same author (not published). In this new approach, similar to Doksum and Hoyland’s work, the failure time is considered to be the rst passage-time of a Gaussian decay process. However, the fatigue processes before and after a change of stress are related without a linear formulation. In this article, we rst describe the model obtained by Bhattacharyya, and discuss the properties of the failure time distribution. Secondly, we present the formulae and the properties of the maximum likelihood estimates for the failure time distribution. Finally, we demonstrate how to apply the step-stress ALT model to elds other than engineering, and illustrate the application of the model studied in this paper to a real data set from a clinical trial.
2. A tampered Brownian motion process model and its properties Let B1 (t) and B2 (t) be two independent Brownian motion processes with positive drifts 1 ¿ 0 and 2 ¿ 0; respectively, and a common di usion parameter 2 . In a
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
17
step-stress ALT environment, a specimen is tested initially under stress v1 and the fatigue grows according to B1 (t). At a preassigned time , if the specimen has not failed yet, the stress is shifted to v2 and the fatigue then grows according to B2 (t) (t¿ ) while starting from the accumulated damage B1 ( ). Denoting B(t) as the fatigue process in the step-stress ALT, the above description can be translated into the following expression: ( t6 ; B1 (t); (1) B(t) = B1 ( ) + B2 (t − ); t ¿ : We call it a tampered Brownian motion process model to re ect the physical meaning of the model: the decay process was interfered with by the switch of the stress. It is assumed that the drift parameters are subject to the relation of 1 ¡ 2 to assure that the fatigue grows faster under the higher stress than the lower one. The failure time of the specimen in the test is then de ned as the rst passage time of the tampered process with respect to a critical boundary , namely, T = inf {t: B(t) ¿ }. Bhattacharyya has derived the probability distribution function for T as 2 1=2 −3=2 ( 2 ) t exp{ −t t6 ; 2 c (1 ; t)}; (2) f(t) = 1=2 −3=2 2 − t ( ) t exp{ 2 [c (1 ; ) + 2 c( ; t)]}s(t); t ¿ ; 2 2
where c(a; b)=1=a−1=b; for a; bp6= 0; 4= c(2 ; 1 ); s(t)=q(t; 4+1; )−q(t; 4−1; ); q(t; a; )=a exp( 12 a2 c( ; t))(a c( ; t)); and (:) is the cumulative distribution function of N(0; 1). The parameters are related to those in the processes by 1 = =1 ; 2 = =2 ; and = 2 =2 : Fig. 1 depicts a path of the tampered Brownian motion process along with the realizations of processes B1 (t) and B2 (t) starting from zero. The corresponding probability density functions of the rst passage time with respect to = 10 are presented in Fig. 2. Let IG(t|; ) denotes an inverse Gaussian distribution with parameters and . Examining the parameters in Eq. (2), we see that 1 and are the parameters in IG(t|1 ; ), and 2 and are those in IG(t|2 ; ) (Chhikara and Folks, 1989). That is, 1 and 2 are the mean lifetime of the specimen under the constant stress v1 and v2 , respectively, and is the common scale parameter for both distributions. It is seen that f(t) on t6 is the pdf of IG(t|1 ; ) restricted on that interval. In contrast, f(t) on t ¿ is di erent from the pdf of IG(t|2 ; ) and involves not only the parameters 2 and but also 1 . This structure re ects the fact that the tampered Brownian motion process on t ¿ is a composition of the increment of B2 (t) on ( ; t] and the cumulated random fatigue of B1 ( ). The distribution function can be shown to have the following features. Proposition 2.1. The probability density function f(t) for the failure time under step-stress ALT is continuous everywhere; and has at least one and at most two local maxima. Namely; f(t) is unimodal if ¡ tm1 and bimodal otherwise; where tm1 =−312 =2+1 (1+912 =42 )1=2 ; the unique mode of IG(t|1 ; ); is the possible mode; and the deÿnite mode is the unique root of the equation =22 + 3=t − 2@ log s(t)=@t = 0:
18
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
Fig. 1. Realizations of the processes. Solid path: a tampered Brownian motion process with 1 = 0:1; 2 = 0:2 and = 0:2. Dotted and dashed paths: Brownian motion processes with = 0:1 and = 0:2, respectively, and a common = 0:2.
Proposition 2.2. There exists ÿnite th moment of the failure time described by the tampered Brownian motion process model for = 1; 2; 3; : : : :
3. Maximum likelihood estimation of the parameters Let = (1 ; 2 ; )0 be the parameters involved in the pdf of (2). Denote f(t) as f(t| ) and s(t) as s(t| ) to clarify the dependence of the function on the parameters. We have ∈ = { ¿ 0; 1 ¿ 0; and 1 ¿ 2 ¿ 0}, which is the parameter space of the distribution family of f(t| ). For n independently observed failure times: t1 ; t2 ; : : : ; tn under step-stress ALT with r (06r6n) of them are observed before time . Without loss of generality, let the rst r ti ’s be in the interval of (0; ]. We have the log-likelihood function: r (t − )2 n P 3P n i 1 log ti − 2 − l( ) = log 2 2 2 i=1 ti 21 i=1 2 1 1 −(n − r) − 2 1 n n P P ti −1 + − 2 log s(ti | ): 22 i=r+1 i=r+1
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
19
Fig. 2. Probability density functions. Solid curve: the pdf of the failure time under step-stress ALT (1 =100, 2 =50 and =2500). Dotted and dashed curves: the pdf’s of the inverse Gaussian distribution IG(t|100; 2500) and IG(t|50; 2500), respectively.
P0 Pn When the extreme case of r = 0 or r = n occurs, the notation i=1 or i=n+1 in the log-likelihood function denotes zero. Expand the quadratic forms and introduce the following statistics: Un =
n P i=1
(ti ∧ );
Vn =
n P i=1
(ti ∧ )−1
and
Wn =
n P
(ti − )+ ;
i=1
where ti ∧ = min{ti ; } and (ti − )+ = max{0; ti − }. The log-likelihood function can be expressed more concisely as n n 3P n log ti + − l( )= log − 2 2 2 i=1 1 2
n P Un Wn + 2 +Vn + log s(ti | ): 2 1 2 i=r+1
(3)
To overcome the diculty caused by the constraint that 1 ¿ 2 , we consider parameter transformation by letting  = k( ), where  = (1 ; 2 ; 3 )0 with 1 =
1 ; 1
2 =
1 1 − ; 2 1
3 = :
(4)
20
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
Di erentiating the transformed log-likelihood function l(Â) with respect to each of the parameters, we have @l(Â) = −Un 1 3 − Wn (1 + 2 )3 + n3 ; @1 n P @s(ti |2 ; 3 )=@2 @l(Â) = −Wn (1 + 2 )3 + ; @2 s(ti |2 ; 3 ) i=r+1 n P @s(ti |2 ; 3 )=@3 n Vn Un 2 Wn @l(Â) 1 − = n1 + − − (1 + 2 )2 + : @3 23 2 2 2 s(ti |2 ; 3 ) i=r+1
From @l(Â)=@1 = 0, 1 can be expressed as the function of 2 by 1 Wn 2 ; 1− (5) 1 = tn n Pn where tn = 1=n i=1 ti . Replacing 1 in the equations of @l(Â)=@2 = 0 and @l(Â)=@3 = 0 by (5), two log-likelihood equations are obtained: W n 3 +
n P Un Wn @s(ti |2 ; 3 )=@2 2 3 − tn =0 n s(ti |2 ; 3 ) i=r+1
ntn +(n−Vn tn )3 −2Wn 2 3 −
n P @s(ti |2 ; 3 )=@3 Un Wn 2 =0: 2 3 +2tn 3 n s(ti |2 ; 3 ) i=r+1
(6)
(7)
In the transformed parameter space 0 , denote the pdf of the failure time with parameter  as f(t|Â) and I (Â) as the Fisher’s information matrix. Let Â0 and 0 be the true parameters in 0 and , respectively. By verifying the conditions and assumptions by Lehmann (1983) for the MLEs in the multiparameter case, we can rst obtain the MLE results for the parameter Â0 (Lehmann, 1983, Theorem 4:1. p. 423) as stated below. Theorem 3.1. For the n iid failure times with t1 ; t2 ; : : : ; tr from (0; ], 06r6n; we have (i) with probability tending to one, there exists a sequence of roots of ˆn of the equations; r @ log f (t |Â) n P P @ log f2 (ti |Â) 1 i + =0 @ @ i=1 i=r+1
(8)
which is consistent for Â0 . √ d (ii) n(ˆn − Â0 ) → N3 (0; I −1 (Â0 )), as n → ∞. (iii) ˆn is asymptotically ecient. The term “n → ∞” in our context always means that the number of observed failures tends to in nity; it does not imply that r is xed or not.
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
21
Note that  = k( ) is a one-to-one and continuously di erentiable transformation, thus Fisher’s information matrix at parameter relates to I (Â) (Lehmann, 1983) by ! −1 0 J ( )I1 (k( ))J 0 ( ) J ( )I2 (k( )) 12 ; where J ( )= 1 I ( )= −1 symmetry I3 (k( )) 2 2 1
2
(9) and Il (k( )) consists of the corresponding components of −E[@2 f(t|Â)=@j @m ]|Â=k( ) . Applying the general delta method to the sequence {ˆn }, which are the MLEs for Â0 with the properties stated in Theorem 3.1, we have the following results. Theorem 3.2. Let ˆ2n and ˆ3n be the roots of the log-likelihood equations (6) and (7) which are calculated from n iid failure times with t1 ; t2 ; : : : ; tr from (0; ], 06r6n. The maximum likelihood estimators ˆ n = (ˆ1n ; ˆ2n ; ˆn )0 for the parameter are given by ˆ1n = The (i) (ii) (iii) (iv)
tn
1 − (Wn =n)ˆ2n
;
ˆ2n =
tn
1 + (tn − Wn =n)ˆ2n
and
ˆn = ˆ3n :
estimators possess the following properties as n → ∞: they are unique with probability tending to one. ˆ n tends to the true parameter 0 strongly. √ ˆ n( n − 0 ) → N3 (0; I −1 ( 0 )) in distribution, where I ( ) is deÿned by (9). √ ˆ n is asymptotically ecient in the sense that n( ˆ jn − j ) → N(0; I −1 ( 0 )j ) in distribution for j = 1; 2; 3.
Following a conventional approach, the covariance matrix of the MLEs ˆ n can be P P−1 estimated by (1=n) ˆ n , where ˆ n = Iˆ( ) = ˆ with Iˆl (k( )) consisting of the corren Pn sponding components of −(1=n) 1 [@2 f(ti |Â)=@j @m ]|Â=k( ˆ ) . n
4. Application of the model Although the partial step-stress ALT model is normally applied in the industrial setting, the underlying structure of the model could be adapted to other contexts, such as the evaluation of drug toxicity where cumulative exposure to the drug will lead to toxicity, presumably at a faster rate with higher doses of drug. If the underlying toxicity is evaluated using a continuous uctuating parameter, such as WBC or SGOT, with respect to a prede ned critical boundary, the tampered Brownian motion process could be used to model time to toxicity data obtained from trials where the dose of the drug is increased in individual subjects over time. To demonstrate our proposal and illustrate the method discussed in this article, we use data taken from a clinical trial of methotrexate (MTX) and trimetrexate in patients with metastatic or recurrent head and neck cancer. On the methotrexate arm patients were to receive a dose of 40 mg=M2 , which could be increased to 60 mg=M2 after
22
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
Fig. 3. White blood cell count (adjusted for baseline) for three patients under administration of MTX. Symbol: (0) no treatment; (1) 40 mg=M2 ; (2) 60 mg=M2 .
day 14 if the WBC remained above 5000. Time to WBC¡ 5000 is the failure time of interest, though for purposes of illustration we will use time to WBC ¡ WBCbaseline − 1000 to accommodate some variation in the baseline value of WBC, which averaged about 6000. Fig. 3 depicts the course of WBC over time for three patients. The rst two patients illustrated had a dose escalation after day 14 and the WBC crossed the critical boundary under the higher dose; the WBC of the third patient dropped below the boundary under the lower dose. The values between measurements are linearly interpolated to obtain the failure time for each individual. A total of 24 patients in the MTX treatment group followed the protocol and failed to maintain normal WBC in the trial. The observed failure times (days) are: 4:67; 9:75; 3:18; 4:67; 12:65; 28:29; 24:5; 5:0; 2:33; 1:94; 16:26; 2:12; 2:92; 10:79; 2:5; 1:79;
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
23
Fig. 4. Histogram of the failure times and the estimated density functions for time to failure under step-stress exposure and under constant exposure to MTX 40 and 60 mg=M2 , respectively.
10:5; 5:8; 3:18; 7:98; 14; 7:48; 2:41; 14. Of these, r =21 failed before =14 days and three failed under the higher dose. The statistics calculated from the data are Un = 171:57, Vn = 5:46, Wn = 27:06, and tn = 8:28: A Newton–Raphson method is applied to solve the equations in Theorem 3.2 with varied initial guesses of 20 = 0:01 (0:01) 0:03 and 30 = 4 (1) 13. For all starting values the roots are ˆ2 = 0:02 and ˆ3 = 9:12, to two decimal places. The MLEs of the three parameters in the tampered Brownian motion process model are estimated as ˆ1 = 8:47; ˆ2 = 7:24; and ˆ = 9:12 following Theorem 3.2, with the estimated variance of 15:68, 19:39 and 0:17, respectively. Fig. 4 presents a histogram of the observed failure times with the estimated pdf of time to failure for the partial step-stress ALT model as well as for the constant stress models. The estimated mean time until failure for the 40 mg=M2 dose is ˆ1 = 8:47 days with a standard error of 3:96, and that for the 60 mg=M2 dose is ˆ2 = 7:24 days with a standard error of 4:40. It is noted that the estimated pdf of the failure time under step-stress ALT is close to that under 40 mg=M2 because the majority of the failures occurred before = 14. In this example, it is not particularly meaningful to speak of the lower dose as the “normal use” condition as compared to the higher dose being a “higher stress” condition. Taken at face value, however, the results from this analysis suggest that the higher dose does not lower the WBC much more rapidly than the lower dose, and could possibly be substituted for the lower dose if greater therapeutic bene t might be derived. The example suggests that the tampered Brownian motion process model could be applicable and useful in some carefully chosen clinical situations.
24
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
Acknowledgements The authors gratefully thank Professor G. K. Bhattacharyya for motivating this research, and are indebted to the editors for their valuable suggestions for simplifying the writing. This work was supported by the National Institute of Health Grant R01-CA45313. References Bhattacharyya, G.K., 1987. Parametric models and inference procedures for accelerated life tests. Invited paper for 46th session of the ISI. Bhattacharyya, G.K., Soejoeti, Z., 1989. A tampered failure rate model for step-stress accelerated life test. Commun. Statist.: Theory Method 18, 1627–1643. Chhikara, R.S., Folks, J.L., 1989. The Inverse Gaussian Distribution — Theory, Methodology, and Applications. Marcel Dekker, Inc, New York. DeGroot, M.H., Goel, P.K., 1979. Bayesian estimation and optimal designs in partially accelerated life testing. Naval Res. Logist. Quart. 26, 223–235. Doksum, K.A., Hoyland, A., 1991. Models for variable-stress accelerated life testing experiments based on Wiener process and inverse Gaussian distribution. Technometrics 34, 74–82. Lehmann, E.L., 1983. Theory of Point Estimation. Wiley, New York. Nelson, W., 1980. Accelerated life testing — step-stress models and data analysis. IEEE Vol.R-29 (2). Schmoyer, R.L., 1991. Nonparametric analysis for two-level single-stress accelerated life tests. Technometrics 33, 175–186. Shaked, M., Singpurwalla, N.D., 1983. Inference for step-stress accelerated life tests. J. Statist. Plann. Inference 7, 295–306.
www.elsevier.com/locate/jspi
A tampered Brownian motion process model for partial step-stress accelerated life testing a Division
Yili Lua; ∗ , Barry Storerb of Statistical and Mathematical Sciences, Lilly Research Laboratories, Eli Lilly and Company, Lilly Corporate Center, Indianapolis, IN 46285, USA b Fred Hutchinson Cancer Research Center, USA
Received 28 July 1998; received in revised form 15 August 2000; accepted 29 August 2000
Abstract In partial step-stress accelerated life testing, models extrapolating data obtained under more severe conditions to infer the lifetime distribution under normal use conditions are needed. Bhattacharyya (Invited paper for 46th session of the ISI, 1987) proposed a tampered Brownian motion process model and later derived the probability distribution from a decay process perspective without linear assumption. In this paper, the model is described and the features of the failure time distribution are discussed. The maximum likelihood estimates of the parameters in the model and their asymptotic properties are presented. An application of models for step-stress accelerated life test to elds other than engineering is described and illustrated by applying the c 2001 Elsevier tampered Brownian motion process model to data taken from a clinical trial. Science B.V. All rights reserved. MSC: 60J65; 62N05 Keywords: Partial acceleration; Brownian motion process; Inverse Gaussian distribution; Maximum likelihood estimation
1. Introduction Accelerated life testing of material is used to collect failure data quickly when the lifetime of a specimen under normal use conditions is too long. The assumption for using an overstress accelerated life test (ALT) is that the lifetime is adversely a ected by some physical conditions which are referred to as stresses. Accelerated life tests which involve higher levels of physical stress than would be experienced under normal conditions will hasten the process of failure. Step-stress ALT is a special kind of ALT in which the stress on each specimen is increased step-by-step over time. When a test ∗
Corresponding author. Tel.: +1 317-277-3954; fax: +1 317-276-5064 E-mail address: lu [email protected] (Y. Lu).
c 2001 Elsevier Science B.V. All rights reserved. 0378-3758/01/$ - see front matter PII: S 0 3 7 8 - 3 7 5 8 ( 0 0 ) 0 0 2 1 5 - 9
16
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
involves two levels of stress with the rst one as the normal one and has a xed time point for changing stress, it is referred to as a partial step-stress ALT. This will be the setting that we deal with throughout this paper, hence, the word partial is omitted hereafter. The goal of statistical inference for step-stress ALT is to estimate the failure behavior of the specimen in conditions of normal use using the failure data obtained under more severe conditions. To achieve this goal, proper statistical models are needed. The pioneering work for step-stress ALT modeling is by DeGroot and Goel (1979) and Nelson (1980). The former proposed a tampered random variable model and the latter a cumulative exposure (CE) model. Alternatively, Bhattacharyya and Soejoeti (1989) suggested a tampered failure rate model, and Doksum and Hoyland (1991) brought forth a model from the decay process point of view. In Doksum and Hoyland’s approach, instead of modeling the failure time distribution directly, a damage growth process is modeled in relation to the acceleration of the stress. The failure time distribution is then obtained as the distribution of the rst passage-time of the process with respect to a critical boundary. In the meantime, nonparametric models have also been proposed (Shaked and Singpurwalla, 1983; Schmoyer, 1991). Most of the parametric models listed above share a common mathematical assumption that there is a linear association between the failure time before and after a change of the stress. For instance, it is assumed that the decay process under the higher stress is linearly associated with the one under the lower stress in Doksum and Hoyland’s model, which yields a tampered random variable model for the failure time with an inverse Gaussian distribution. It can also be shown that the tampered random variable models compose a proper subset of the cumulative exposure (CE) models, when the CE models are limited to those within location-scale families. Therefore, the availability of non-linear models is very limited when the linear assumption fails in the real world. This paper examines another approach for modeling step-stress ALT. The idea of the approach was mentioned by Bhattacharyya (1987) and the model was later derived by the same author (not published). In this new approach, similar to Doksum and Hoyland’s work, the failure time is considered to be the rst passage-time of a Gaussian decay process. However, the fatigue processes before and after a change of stress are related without a linear formulation. In this article, we rst describe the model obtained by Bhattacharyya, and discuss the properties of the failure time distribution. Secondly, we present the formulae and the properties of the maximum likelihood estimates for the failure time distribution. Finally, we demonstrate how to apply the step-stress ALT model to elds other than engineering, and illustrate the application of the model studied in this paper to a real data set from a clinical trial.
2. A tampered Brownian motion process model and its properties Let B1 (t) and B2 (t) be two independent Brownian motion processes with positive drifts 1 ¿ 0 and 2 ¿ 0; respectively, and a common di usion parameter 2 . In a
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
17
step-stress ALT environment, a specimen is tested initially under stress v1 and the fatigue grows according to B1 (t). At a preassigned time , if the specimen has not failed yet, the stress is shifted to v2 and the fatigue then grows according to B2 (t) (t¿ ) while starting from the accumulated damage B1 ( ). Denoting B(t) as the fatigue process in the step-stress ALT, the above description can be translated into the following expression: ( t6 ; B1 (t); (1) B(t) = B1 ( ) + B2 (t − ); t ¿ : We call it a tampered Brownian motion process model to re ect the physical meaning of the model: the decay process was interfered with by the switch of the stress. It is assumed that the drift parameters are subject to the relation of 1 ¡ 2 to assure that the fatigue grows faster under the higher stress than the lower one. The failure time of the specimen in the test is then de ned as the rst passage time of the tampered process with respect to a critical boundary , namely, T = inf {t: B(t) ¿ }. Bhattacharyya has derived the probability distribution function for T as 2 1=2 −3=2 ( 2 ) t exp{ −t t6 ; 2 c (1 ; t)}; (2) f(t) = 1=2 −3=2 2 − t ( ) t exp{ 2 [c (1 ; ) + 2 c( ; t)]}s(t); t ¿ ; 2 2
where c(a; b)=1=a−1=b; for a; bp6= 0; 4= c(2 ; 1 ); s(t)=q(t; 4+1; )−q(t; 4−1; ); q(t; a; )=a exp( 12 a2 c( ; t))(a c( ; t)); and (:) is the cumulative distribution function of N(0; 1). The parameters are related to those in the processes by 1 = =1 ; 2 = =2 ; and = 2 =2 : Fig. 1 depicts a path of the tampered Brownian motion process along with the realizations of processes B1 (t) and B2 (t) starting from zero. The corresponding probability density functions of the rst passage time with respect to = 10 are presented in Fig. 2. Let IG(t|; ) denotes an inverse Gaussian distribution with parameters and . Examining the parameters in Eq. (2), we see that 1 and are the parameters in IG(t|1 ; ), and 2 and are those in IG(t|2 ; ) (Chhikara and Folks, 1989). That is, 1 and 2 are the mean lifetime of the specimen under the constant stress v1 and v2 , respectively, and is the common scale parameter for both distributions. It is seen that f(t) on t6 is the pdf of IG(t|1 ; ) restricted on that interval. In contrast, f(t) on t ¿ is di erent from the pdf of IG(t|2 ; ) and involves not only the parameters 2 and but also 1 . This structure re ects the fact that the tampered Brownian motion process on t ¿ is a composition of the increment of B2 (t) on ( ; t] and the cumulated random fatigue of B1 ( ). The distribution function can be shown to have the following features. Proposition 2.1. The probability density function f(t) for the failure time under step-stress ALT is continuous everywhere; and has at least one and at most two local maxima. Namely; f(t) is unimodal if ¡ tm1 and bimodal otherwise; where tm1 =−312 =2+1 (1+912 =42 )1=2 ; the unique mode of IG(t|1 ; ); is the possible mode; and the deÿnite mode is the unique root of the equation =22 + 3=t − 2@ log s(t)=@t = 0:
18
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
Fig. 1. Realizations of the processes. Solid path: a tampered Brownian motion process with 1 = 0:1; 2 = 0:2 and = 0:2. Dotted and dashed paths: Brownian motion processes with = 0:1 and = 0:2, respectively, and a common = 0:2.
Proposition 2.2. There exists ÿnite th moment of the failure time described by the tampered Brownian motion process model for = 1; 2; 3; : : : :
3. Maximum likelihood estimation of the parameters Let = (1 ; 2 ; )0 be the parameters involved in the pdf of (2). Denote f(t) as f(t| ) and s(t) as s(t| ) to clarify the dependence of the function on the parameters. We have ∈ = { ¿ 0; 1 ¿ 0; and 1 ¿ 2 ¿ 0}, which is the parameter space of the distribution family of f(t| ). For n independently observed failure times: t1 ; t2 ; : : : ; tn under step-stress ALT with r (06r6n) of them are observed before time . Without loss of generality, let the rst r ti ’s be in the interval of (0; ]. We have the log-likelihood function: r (t − )2 n P 3P n i 1 log ti − 2 − l( ) = log 2 2 2 i=1 ti 21 i=1 2 1 1 −(n − r) − 2 1 n n P P ti −1 + − 2 log s(ti | ): 22 i=r+1 i=r+1
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
19
Fig. 2. Probability density functions. Solid curve: the pdf of the failure time under step-stress ALT (1 =100, 2 =50 and =2500). Dotted and dashed curves: the pdf’s of the inverse Gaussian distribution IG(t|100; 2500) and IG(t|50; 2500), respectively.
P0 Pn When the extreme case of r = 0 or r = n occurs, the notation i=1 or i=n+1 in the log-likelihood function denotes zero. Expand the quadratic forms and introduce the following statistics: Un =
n P i=1
(ti ∧ );
Vn =
n P i=1
(ti ∧ )−1
and
Wn =
n P
(ti − )+ ;
i=1
where ti ∧ = min{ti ; } and (ti − )+ = max{0; ti − }. The log-likelihood function can be expressed more concisely as n n 3P n log ti + − l( )= log − 2 2 2 i=1 1 2
n P Un Wn + 2 +Vn + log s(ti | ): 2 1 2 i=r+1
(3)
To overcome the diculty caused by the constraint that 1 ¿ 2 , we consider parameter transformation by letting  = k( ), where  = (1 ; 2 ; 3 )0 with 1 =
1 ; 1
2 =
1 1 − ; 2 1
3 = :
(4)
20
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
Di erentiating the transformed log-likelihood function l(Â) with respect to each of the parameters, we have @l(Â) = −Un 1 3 − Wn (1 + 2 )3 + n3 ; @1 n P @s(ti |2 ; 3 )=@2 @l(Â) = −Wn (1 + 2 )3 + ; @2 s(ti |2 ; 3 ) i=r+1 n P @s(ti |2 ; 3 )=@3 n Vn Un 2 Wn @l(Â) 1 − = n1 + − − (1 + 2 )2 + : @3 23 2 2 2 s(ti |2 ; 3 ) i=r+1
From @l(Â)=@1 = 0, 1 can be expressed as the function of 2 by 1 Wn 2 ; 1− (5) 1 = tn n Pn where tn = 1=n i=1 ti . Replacing 1 in the equations of @l(Â)=@2 = 0 and @l(Â)=@3 = 0 by (5), two log-likelihood equations are obtained: W n 3 +
n P Un Wn @s(ti |2 ; 3 )=@2 2 3 − tn =0 n s(ti |2 ; 3 ) i=r+1
ntn +(n−Vn tn )3 −2Wn 2 3 −
n P @s(ti |2 ; 3 )=@3 Un Wn 2 =0: 2 3 +2tn 3 n s(ti |2 ; 3 ) i=r+1
(6)
(7)
In the transformed parameter space 0 , denote the pdf of the failure time with parameter  as f(t|Â) and I (Â) as the Fisher’s information matrix. Let Â0 and 0 be the true parameters in 0 and , respectively. By verifying the conditions and assumptions by Lehmann (1983) for the MLEs in the multiparameter case, we can rst obtain the MLE results for the parameter Â0 (Lehmann, 1983, Theorem 4:1. p. 423) as stated below. Theorem 3.1. For the n iid failure times with t1 ; t2 ; : : : ; tr from (0; ], 06r6n; we have (i) with probability tending to one, there exists a sequence of roots of ˆn of the equations; r @ log f (t |Â) n P P @ log f2 (ti |Â) 1 i + =0 @ @ i=1 i=r+1
(8)
which is consistent for Â0 . √ d (ii) n(ˆn − Â0 ) → N3 (0; I −1 (Â0 )), as n → ∞. (iii) ˆn is asymptotically ecient. The term “n → ∞” in our context always means that the number of observed failures tends to in nity; it does not imply that r is xed or not.
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
21
Note that  = k( ) is a one-to-one and continuously di erentiable transformation, thus Fisher’s information matrix at parameter relates to I (Â) (Lehmann, 1983) by ! −1 0 J ( )I1 (k( ))J 0 ( ) J ( )I2 (k( )) 12 ; where J ( )= 1 I ( )= −1 symmetry I3 (k( )) 2 2 1
2
(9) and Il (k( )) consists of the corresponding components of −E[@2 f(t|Â)=@j @m ]|Â=k( ) . Applying the general delta method to the sequence {ˆn }, which are the MLEs for Â0 with the properties stated in Theorem 3.1, we have the following results. Theorem 3.2. Let ˆ2n and ˆ3n be the roots of the log-likelihood equations (6) and (7) which are calculated from n iid failure times with t1 ; t2 ; : : : ; tr from (0; ], 06r6n. The maximum likelihood estimators ˆ n = (ˆ1n ; ˆ2n ; ˆn )0 for the parameter are given by ˆ1n = The (i) (ii) (iii) (iv)
tn
1 − (Wn =n)ˆ2n
;
ˆ2n =
tn
1 + (tn − Wn =n)ˆ2n
and
ˆn = ˆ3n :
estimators possess the following properties as n → ∞: they are unique with probability tending to one. ˆ n tends to the true parameter 0 strongly. √ ˆ n( n − 0 ) → N3 (0; I −1 ( 0 )) in distribution, where I ( ) is deÿned by (9). √ ˆ n is asymptotically ecient in the sense that n( ˆ jn − j ) → N(0; I −1 ( 0 )j ) in distribution for j = 1; 2; 3.
Following a conventional approach, the covariance matrix of the MLEs ˆ n can be P P−1 estimated by (1=n) ˆ n , where ˆ n = Iˆ( ) = ˆ with Iˆl (k( )) consisting of the corren Pn sponding components of −(1=n) 1 [@2 f(ti |Â)=@j @m ]|Â=k( ˆ ) . n
4. Application of the model Although the partial step-stress ALT model is normally applied in the industrial setting, the underlying structure of the model could be adapted to other contexts, such as the evaluation of drug toxicity where cumulative exposure to the drug will lead to toxicity, presumably at a faster rate with higher doses of drug. If the underlying toxicity is evaluated using a continuous uctuating parameter, such as WBC or SGOT, with respect to a prede ned critical boundary, the tampered Brownian motion process could be used to model time to toxicity data obtained from trials where the dose of the drug is increased in individual subjects over time. To demonstrate our proposal and illustrate the method discussed in this article, we use data taken from a clinical trial of methotrexate (MTX) and trimetrexate in patients with metastatic or recurrent head and neck cancer. On the methotrexate arm patients were to receive a dose of 40 mg=M2 , which could be increased to 60 mg=M2 after
22
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
Fig. 3. White blood cell count (adjusted for baseline) for three patients under administration of MTX. Symbol: (0) no treatment; (1) 40 mg=M2 ; (2) 60 mg=M2 .
day 14 if the WBC remained above 5000. Time to WBC¡ 5000 is the failure time of interest, though for purposes of illustration we will use time to WBC ¡ WBCbaseline − 1000 to accommodate some variation in the baseline value of WBC, which averaged about 6000. Fig. 3 depicts the course of WBC over time for three patients. The rst two patients illustrated had a dose escalation after day 14 and the WBC crossed the critical boundary under the higher dose; the WBC of the third patient dropped below the boundary under the lower dose. The values between measurements are linearly interpolated to obtain the failure time for each individual. A total of 24 patients in the MTX treatment group followed the protocol and failed to maintain normal WBC in the trial. The observed failure times (days) are: 4:67; 9:75; 3:18; 4:67; 12:65; 28:29; 24:5; 5:0; 2:33; 1:94; 16:26; 2:12; 2:92; 10:79; 2:5; 1:79;
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
23
Fig. 4. Histogram of the failure times and the estimated density functions for time to failure under step-stress exposure and under constant exposure to MTX 40 and 60 mg=M2 , respectively.
10:5; 5:8; 3:18; 7:98; 14; 7:48; 2:41; 14. Of these, r =21 failed before =14 days and three failed under the higher dose. The statistics calculated from the data are Un = 171:57, Vn = 5:46, Wn = 27:06, and tn = 8:28: A Newton–Raphson method is applied to solve the equations in Theorem 3.2 with varied initial guesses of 20 = 0:01 (0:01) 0:03 and 30 = 4 (1) 13. For all starting values the roots are ˆ2 = 0:02 and ˆ3 = 9:12, to two decimal places. The MLEs of the three parameters in the tampered Brownian motion process model are estimated as ˆ1 = 8:47; ˆ2 = 7:24; and ˆ = 9:12 following Theorem 3.2, with the estimated variance of 15:68, 19:39 and 0:17, respectively. Fig. 4 presents a histogram of the observed failure times with the estimated pdf of time to failure for the partial step-stress ALT model as well as for the constant stress models. The estimated mean time until failure for the 40 mg=M2 dose is ˆ1 = 8:47 days with a standard error of 3:96, and that for the 60 mg=M2 dose is ˆ2 = 7:24 days with a standard error of 4:40. It is noted that the estimated pdf of the failure time under step-stress ALT is close to that under 40 mg=M2 because the majority of the failures occurred before = 14. In this example, it is not particularly meaningful to speak of the lower dose as the “normal use” condition as compared to the higher dose being a “higher stress” condition. Taken at face value, however, the results from this analysis suggest that the higher dose does not lower the WBC much more rapidly than the lower dose, and could possibly be substituted for the lower dose if greater therapeutic bene t might be derived. The example suggests that the tampered Brownian motion process model could be applicable and useful in some carefully chosen clinical situations.
24
Y. Lu, B. Storer / Journal of Statistical Planning and Inference 94 (2001) 15–24
Acknowledgements The authors gratefully thank Professor G. K. Bhattacharyya for motivating this research, and are indebted to the editors for their valuable suggestions for simplifying the writing. This work was supported by the National Institute of Health Grant R01-CA45313. References Bhattacharyya, G.K., 1987. Parametric models and inference procedures for accelerated life tests. Invited paper for 46th session of the ISI. Bhattacharyya, G.K., Soejoeti, Z., 1989. A tampered failure rate model for step-stress accelerated life test. Commun. Statist.: Theory Method 18, 1627–1643. Chhikara, R.S., Folks, J.L., 1989. The Inverse Gaussian Distribution — Theory, Methodology, and Applications. Marcel Dekker, Inc, New York. DeGroot, M.H., Goel, P.K., 1979. Bayesian estimation and optimal designs in partially accelerated life testing. Naval Res. Logist. Quart. 26, 223–235. Doksum, K.A., Hoyland, A., 1991. Models for variable-stress accelerated life testing experiments based on Wiener process and inverse Gaussian distribution. Technometrics 34, 74–82. Lehmann, E.L., 1983. Theory of Point Estimation. Wiley, New York. Nelson, W., 1980. Accelerated life testing — step-stress models and data analysis. IEEE Vol.R-29 (2). Schmoyer, R.L., 1991. Nonparametric analysis for two-level single-stress accelerated life tests. Technometrics 33, 175–186. Shaked, M., Singpurwalla, N.D., 1983. Inference for step-stress accelerated life tests. J. Statist. Plann. Inference 7, 295–306.