Extrapolation and confounding in accelerated material degradation and failure studies

Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

Contents lists available at ScienceDirect

Journal of Statistical Planning and Inference journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / j s p i

Extrapolation and confounding in accelerated material degradation and failure studies夡 M.J. LuValle OFS Laboratories, USA

A R T I C L E

I N F O

Available online 28 May 2008 Keywords: Extrapolation Kinetics Confounding Risk orthogonal experiments Sequential experiments

A B S T R A C T

The fundamental difficulty with inference in nontrivial extrapolation where model selection is involved from a rich space of models is that any model estimated in one regime used for decision making in another is fundamentally confounded with disruptive alternatives. These are alternative models which if true would support a diametrically opposed action from the one the estimated model supports. One strategy to support extrapolation and reduce arbitrary fitting and confounding is to force the model to derive from the same mathematical structure that underlies the substantive science appropriate for the phenomena. Then statistical model fitting follows the form of theory generation in artificial intelligence, with statistical model selection tools and the statistician taking the place of the inference engine. However, the problem of confounding still remains. Activity mapping is a recent development which, for any given specific form of alternative kinetic model, identifies the confounded models having that form. In this paper, it is shown how activity maps may be used to help design experiments with improved power for detecting confounding disruptive alternatives in the context of fitting kinetic models to material degradation or failure data with the goal of extrapolation. © 2008 Elsevier B.V. All rights reserved.

1. Introduction 1.1. Overview Usually the goal of accelerated testing for new technologies is to obtain an estimate of the field reliability of a new device or system without having to wait for field reliability data. For the purpose of this paper we confine our attention to accelerated testing of material systems, with concern only with physical failure and degradation modes. The introduction of a significantly new material system, or application of an existing material system in a new environment, can pose significant risk. Even with accelerated testing, when field experience is lacking, failure modes/mechanisms can elude detection. Alternatively, failure modes observed in accelerated tests may in fact be artificial, but not reveal this aspect of their nature at highly accelerated test levels. Stress voiding in aluminum lines on silicon integrated circuits (Turner and Wendel, 1985; Yost et al., 1989), and a secondary relaxation mode in optically induced gratings in boron doped optical fibers (LuValle et al., 1998) are examples of the first. The conductive anodic filament (CAF) failure mode (Welsher et al., 1980; Augis et al., 1989; LuValle and Hines, 1992; Meeker and LuValle, 1995) is an example of the second. Recent theoretical results (LuValle, 2007) provide a framework for identifying processes that provide a sort of extreme confounding with an existing model based on accelerated test data. The confounding processes are called evanescent processes (LuValle et al., 2002, 2004; LuValle, 2004a, b), because they vanish, or evanesce at highly accelerated conditions, but dominate at 夡 Work partially supported by funds from the Naval Research Laboratory. E-mail address: [email protected] 0378-3758/$ - see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2008.05.035

1708

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

operating conditions. The purpose of this paper is to expand on the work that allowed identification of the confounding models (LuValle, 2007) to develop some statistical theory of experiment design, allowing tradeoffs of design vs. approximate estimated power of detecting the confounding models to be made. The remainder of this section provides a brief introduction to kinetic modeling, evanescent processes, and activity maps. In the next section we consider two situations, one where it is possible to design a special kind of experiment (LuValle et al., 2002, 2004; LuValle, 2004a), called risk orthogonal, where the test statistic does not require consideration of extrapolation variance, and the other where a sequential approach to design, trading off variance and bias in the test statistic is required. The variance reduction component of the latter derives from the same considerations (Elfving, 1952; Chermoff, 1972) that others have used to derive optimal and robust designs for accelerated testing under the assumption that the form of the accelerated testing model is known (e.g., Meeker and Nelson, 1975; Meeker and Hahn 1985; Nelson, 1990). A simple example of a risk orthogonal experiment is given at the end of the second section. In the third section, a numerical study using simple approximations is made of a particular case where risk orthogonal alternatives are not available (LuValle et al., 2006). The final section summarizes the results and reviews open questions about the approach.

1.2. Kinetic modeling and evanescent processes There are many books providing an overview of accelerated testing, either as chapters, or as the main purpose of the book (e.g., Viertl, 1988; Nelson, 1990; Meeker and Escobar, 1998; Bagdonavicius and Nikulin, 2002). While some of the models are tied to physical models, it is not a primary tenet for a model to be discussed. The accelerated life model is the canonical model that ties to a physical notion (LuValle et al., 1988, 2004; Klinger, 1992), but to the author's knowledge generalization of that to arbitrary kinetic models has been followed only in one, more specialized, book (LuValle et al., 2004), and in few publications (LuValle, 1999, 2004a). Some misapprehensions exist about kinetic modeling. Usually it is thought that it is necessary to have highly specialized knowledge about failure modes before kinetic modeling can be applied. The author's approach is to consider the mathematical structure provided by kinetic models as a modeling language, which can be used to generate reasonable hypotheses about some of the structure that data will follow. This is motivated in part by some work in artificial intelligence on knowledge discovery (Lenat and Brown, 1984). Essentially, the statistician serves as an inference engine in the syntactic framework setup by the software, and the product is a model of the degradation or failure process, consistent with using a kinetic analogy to describe all degradation and failure processes. For many degradation and failure processes, the kinetics framework is exactly correct. A book and corresponding freeware are focused on this approach (LuValle et al., 2004). The difficulty in this approach lies in the fact that kinetic models deal with changes in a (usually) unobserved state vector, so it is necessary to assume some form of translation function leading from that state vector to the observable. For the accelerated life model with failure time data, the default assumption is that each device possesses a random threshold, one that does not change with stress, so that if some function of the state vector exceeds a threshold failure occurs. This quite simply generalizes to any other kinetic model. In some cases (e.g. electromigration and stress voiding in integrated circuits, and radiation and hydrogen induced loss in optical fibers) the translation function is relatively well defined by device physics. In other cases, simple assumptions may serve to start. Experimental or theoretical verification of such assumptions is a task that must be coordinated with the experimentalist/subject matter specialist. Kinetic models can give quite interesting alternatives to the standard model, and naturally provide solutions to the problem of understanding step or changing stress experiments. For example, the limited failure probability model (Meeker, 1987) in which only a portion of the devices will ever fail arises naturally in many situations, and in some the proportion (ever) failing changes with changing stress. However, the original formulation of the model results in a stable proportion failing. There are two very simple kinetic models which result in a limited failure probability model, both of which naturally scale the proportion failing with changing stress (LuValle and Hines, 1992; LuValle et al., 2004). One, a branching process results in stress hardening, whereas the other, a reversible process does not. Another example, when dealing with a random material such as a glass, many processes can have a distribution of activation energies across sites in the glass (Primak, 1955). This results in a power law relation in absolute temperature for scaling between stresses (LuValle et al., 1998, 2004). More specifically, it is useful to understand what the space of kinetic models looks like. One description given before (LuValle, 2004b) has the form:

jy1 = f1 (y1 , . . . , yk ) jt ...

jyk = fk (y1 , . . . , yk ), jt

(1)

where the fi take the form of polynomials in the y1 , . . . , yk . The yi typically represent concentrations in standard chemical kinetic texts, but can be more general (LuValle et al., 1988).

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

1709

In the general case the fi take the form: fi =

ri 

⎛ ⎞ k  zlij Cij ⎝ y ⎠ . l

j=1

(2)

l=1

The expression on the right-hand side of (2) is a polynomial. The Cij are the rate constants for each elementary process resulting in the increment or decrement of yi . Thus for thermally accelerated reactions the Cij will be Arrhenius terms. The term ri counts the number of elementary processes which effect the concentration of the yi . The vector zij = (z1ij , . . . , zkij ) lies on the k-fold cross product of the positive integers (Z + )k . An alternative representation is to separate (1) and (2) into a linear portion (first and zeroth-order kinetics) where the zij = (z1ij , . . . , zkij ) have the form of the canonical basis vectors (k–1 zeros and 1 one), and a nonlinear portion. A simple one-step model corresponding to the accelerated life model can have the form: dA1 = − k1 An1 , dt dA2 = k1 An1 . dt

(3)

Here the Ai represent the concentration of species i, different species, failure corresponds to the accumulation of A2 beyond a random (fixed by the experimental unit) threshold, and n is termed the order of the reaction. The “rate constants” (e.g. k1 ) for the systems of differential equations depend on the stress in certain reasonably constrained ways. For example, the Arrhenius model is often a good approximation for the dependence of an elementary kinetic process on temperature. The Arrhenius model has the form  Ea , ki =  exp − kT

(4)

where Ea is called the activation energy, and represents a barrier to the reaction, k is Boltzmanns constant and T is absolute temperature. Typically  is assumed constant across temperature. However, investigations into the theory of the rates of chemical reactions (Glasstone et al., 1941; Laidler, 1987) show that at least in simpler cases  ∝ (kT/ h) ¯ exp((S)) where h¯ is Planck's constant, and (S) is the difference in entropy between the activated state, and the reactants. In addition to the obvious linear relation, there is also a temperature dependence embedded in the entropy that is dependent on the details of the reaction. However, for a reasonable temperature range (e.g. 0--500 ◦ C) (4) is sufficiently accurate. The proportionality constant is particular to the reaction, but independent of temperature (as in e.g. Matthewson, 1998). An equilibration based model which results in the limited failure probability model has the form: dA1 = − k 1 A1 + k 2 A2 , dt dA2 = k 1 A1 − k 2 A2 , dt

(5)

again where failure corresponds to the accumulation of A2 beyond a random (fixed by the experimental unit) threshold. A branching model with exactly the same integrated law for constant stress experiments has the form: dA1 = − (k1 + k2 )A1 , dt dA2 = k 1 A1 , dt dA3 = k 2 A1 . dt

(6)

The difference between the two models lies in the fact that under fluctuating stress, (6) as an evanescent process can result in stress hardening of individual devices, while (5) will not. Even though there are a large number of evanescent processes possible for any postulated base kinetic process (LuValle, 2004b), the processes given by (3), (5), and (6) are sufficiently complex to produce disruptive evanescent processes (processes, that once known about, reverse the conclusion of the extrapolation, but which are undetectable in the original data) for any base kinetic process in the following sense: • If the original extrapolation indicates the new material system is reliable, then any of the processes given by systems (3), (5), or (6) added as entirely independent processes can specify a disruptive evanescent process. • If the original extrapolation indicates the new material system is unreliable then (5) or (6), integrated into the original process are sufficient to produce disruptive evanescent processes.

1710

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

20

1 0.8 0.6 0.4 0.2 0 20

10

15

2

10

5

log10 (n

0 u)

-5

-10

0.5

1

2.5

3

log10 (nu)

Z

15

1.5 a E

5

0

-5

-10 0.0

0.5

1.0

1.5 Ea

2.0

2.5

3.0

Fig. 1. Integrated activity for 25 years at 25 ◦ C, calculated from (7).

20

15

10 0.8 0.6 0.4 0.2 0 20

log10 (nu)

Z

1

5

15

10

5 log10

0 (nu)

-5

-10

0.5

1

1.5

2

2.

5

3 0

Ea -5

-10 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Ea Fig. 2. Integrated activity for 6 h, 300 ◦ C calculated from (7).

1.3. Activity mapping Given a form for a particular evanescent process, say for example Eq. (3), with the rate constant depending on temperature through the Arrhenius relationship, then the activity map is a map of the maximum normalized integrated chemical activity of the process (normalized with respect to the maximum integrated activity possible for hypothesized species of concern) over the

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

1711

20

15

10

log10 (nu)

1

Z

0.5 0 -0.5

5

-1

20

15

5

1.

10

1

5

log1

0

0 (nu

)

-5

2

2.

Ea

5

3 0

5

0.

-10

-5

-10 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Ea Fig. 3. Activity map, subtracting Fig. 1 from Fig. 2.

course of the experiment minus the maximum normalized integrated chemical activity of the process during life. Another term for maximum normalized integrated chemical activity, in the special case where no reaction has occurred prior to the time of interest, is extent of reaction. In this paper, the term integrated activity, or activity will be used interchangeably with maximum normalized integrated chemical activity. In the case of (3), with an Arrhenius dependence on temperature with n = 1, for life and experiments without thermal pretreatment (e.g. burn in) the normalized integrated solution is just A2t A10 = A10 A10

   Ea . 1 − exp −t exp − kT

(7)

Evaluated for each (, Ea ) pair for a given time and absolute temperature of the experiment. For a life of 25 years, at a temperature of 25 ◦ C, this integrated activity takes the form shown in Fig. 1. The perspective plot on the left side of Fig. 1 shows (7), the integrated activity during life plotted on the z axis with log10 () in Hertz on the y-axis and Ea in electron volts on the x-axis. The plot on the right side is an image plot showing the different levels of z. On these axis, (7) is very sharp. For an experiment of 6 h at 300 ◦ C the integrated activity takes the form shown in Fig. 2. The activity map (Fig. 2 minus Fig. 1) is shown in Fig. 3. The negative portion of the activity map (below 0 on the perspective plot, or the light region on the image map) is the region where the activity during life is higher than that during the experiment. So this is the region where potential disruptive alternatives having the form given in (3) with n = 1 may exist. In particular, given that the range of the parameters of the process are as large as necessary (which can be justified by physical argument), and making the approximation that any failure mode in regions of overlap would have given an unmistakeable signal (which can be quantified by statistical argument based on how representative the sample used in the base accelerated testing is) the volume of this region (hyper volume in more complex processes) is precisely the appropriate support for a prior, for a disruptive alternative of this form. Within this paper, that region will correspond to the volume of uncertainty. Note that this region is determined only by the form of the experiments run to date, and does not require observed degradation or failure to date, exactly as one would expect for studying confounding. To see how the volume of uncertainty is affected by further experiments, we consider Fig. 4, where we see the effect of both the experiment of 6 h at 300 ◦ C, and an additional experiment of 1 year at 100 ◦ C. Focus on the negative portion of the perspective plot and the lighter part of the image plot. This negative region is where integrated activity during life exceeds integrated activity during the experiments. Both the width of the volume of uncertainty across the parameter space, and the maximum depth are significantly reduced from Fig. 3. To see the reduction in depth note that the negative part of the perspective plot in Fig. 3

1712

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

20

15

10

log10 (nu)

1

Z

0.5 0 -0.5 -1 20

5

15

5

10 log10

5 (nu)

1 0

-5

1.

2

2.

5

3 0

Ea

5

0. -10

-5

-10 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Ea Fig. 4. Activity map of life vs. combined experiments 6 h at 300 ◦ C and 1 year at 100 ◦ C. Reduction of volume of uncertainty by 71% from Fig. 3.

reaches down to −1 on the z-axis, where as the negative part of the perspective plot in Fig. 4 terminates well above −1. Similarly, the width of the light triangle of the image plot in Fig. 3 is significantly broader than the width of the light triangle in Fig. 4. The reduction in volume was ∼71%. The map in Fig. 4 was constructed by taking the maximum of the integrated activity between the two experiments, and subtracting the integrated activity during life at each point in the parameter space. The activity map is a primitive version of an approximate tool, the demarcation map (LuValle et al., 1998). In fact, the demarcation map approximates the boundaries of the activity map. This relationship is discussed more closely in LuValle (2007). While the approximation is quite good for simple cases, it becomes very difficult to calculate as the complexity of the process increases. Thus, while model (3) gives a simple analytical calculation for the demarcation approximation, model (5) results in a very complex numerical calculation to even define a demarcation approximation. In contrast, it is quite straightforward to set up the calculations for the activity map for the model given by (5). To make use of the information provided by the activity map to design experiments, it is necessary to understand the operating characteristics of particular patterns of design coupled with inference. We will discuss two situations in the next section. One where it is possible to use the activity map to define “risk orthogonal” experiments, for which the inference problem is simply to detect some new signal against an essentially zero signal from the base model, and the other where a sequential approach is appropriate for improving the power of detecting a departure from a non-zero estimated model. An example of identifying risk orthogonal experimental conditions will be given at the end of the second section. In the third section a numerical study of a real problem using the sequential approach is presented. 2. Risk orthogonal experiments vs. sequential experiments 2.1. Conceptual discussion: risk orthogonal vs. sequential experiments As mentioned above, from the statistical viewpoint we will consider two cases, but before proceeding to that, we discuss in what situations those cases arise. It is possible to classify the decision problems into four situations listed below: I. No degradation or failure is observed during the accelerated testing. In this case extrapolation indicates the product is reliable, and any further experiments will be done against testing a null model of no failure or degradation, thus this situation already results in a situation where the experiments are risk orthogonal, without special work.

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

1713

II. Degradation or failure is observed, but extrapolation clearly indicates that the product is reliable. In this case, the disruptive alternative would be that there is confounded model which results in significantly lower reliability. Here depending on the natural structural model for the phenomena in question, it may be possible to use activity to define risk orthogonal experiments. In all cases, in this situation, the hypothesized evanescent process can be designed to be structurally independent (LuValle, 2007) (no connection between the states involved in the evanescent processes and the states involved in the initial degradation process). III. Extrapolation clearly indicates that the product is unreliable. In this case the disruptive alternative would be that there is a confounded model which results in significantly higher reliability. To see how such can happen consider the systems (5) and (6) from the previous section. If k1 ?k2 in the accelerated stress regime, but k1 >k2 in operating regime, it can happen that in the accelerated regime, it is not apparent that there is a limited failure probability model, but in the operating regime, no failures occur. One example of this would be if all tests were at very high temperature, k1 had a very high activation energy, and k2 had a very low activation energy, with respective values of the frequency factor ( from (4)) that allow a trade off. Investigation in this case usually only occurs when there is enormous economic loss in abandoning the technology direction. This was the case with the conductive anodic filament failure mode (Welsher et al., 1980; Augis et al., 1989; LuValle and Hines, 1992; Meeker and LuValle, 1995). In this case, it is inherently impossible to construct risk orthogonal design. However, some times analysis tools and experimental designs can seize on characteristics of the alternative model that can be tested in new experiments independently of the extrapolation from the first model. (E.g. testing for curvature in a quantile--quantile plot, where we plot quantiles from one stress condition against quantiles corresponding to equivalent probability of failure from another stress condition.) IV. Extrapolation that does not provide a clear indication of the reliability of the product. Then the initial effort needs to be improving precision of the extrapolation of reliability, using insight into design and analysis gathered from the second and third situation. Although risk orthogonal experiments have been defined elsewhere (LuValle et al., 2004; LuValle, 2004a, b), it was done before the discovery of activity maps. With activity mapping, it is possible to give a precise operating definition of a risk orthogonal experiment in terms of confidence sets for the parameters of the process, and the value of an activity map for comparing the initial experiments and final experiments. A new experiment set {E1∗ , . . . , Ek∗ } is defined to be (,) risk orthogonal with respect to a failure or degradation mechanism

modeled using parameter set , based on experiments {E1 , . . . , En } if the maximum integrated activity over the new experiment set {E1∗ , . . . , Ek∗ } of any process within a specified 1 −  confidence region of  is less than . With any estimated model, there will be an assignment of integrated activity to the observations of the prior experiments, either as thresholds for failure time, or as values assigned to observed degradation. Thus,  can be translated to either a probability of a failure, or the amount of degradation, based on the original modeling. In the case that some degree of risk orthogonality is possible for new experiments probing the alternatives in question, statistical tests and power comparisons for experiment designs are greatly simplified. To be specific, suppose that X1 , . . . , Xn are the relevant covariates of the experiments to date, and f (, Xi ) is the relevant integrated kinetic model, including the translation from the state vector to the observable. In particular, f (, X) = g(, A( , X)), where A is the integrated activity, g is the translation function and  = (, ) the concatenation of the two parameter vectors. Assume g has bounded derivatives over the range [0,1]. Assume a sequence {Xj } corresponding to the further experiments. Consider a test at point j. Assuming the estimated j=1, ...,k 2 ˆ error variance is , with additive homogeneous error, the test statistic at a given point depends asymptotically on an error variance like: ⎛ ⎞ ⎛



 ⎞−1 ˆ , X ) ˆ, X ) ˆ ˆ j f (  j f (  j f (  , X ) j f (  , X ) j ⎝ j ⎟ i i ⎠

ˆ 2 ⎜ + 1⎠ (8) ⎝ j j i=1, ...,n j i=1, ...,n j to test the difference between the extrapolation and the measured data. Here V  is the transpose of V. Assume additive Gaussian measurement error. By definition in the (, ) risk orthogonal case, we know with confidence 1 − , the value of A( , Xj ) < . The bounded derivative of g implies that g satisfies the Lipschitz condition. Thus if we wish to test the hypothesis that the model currently selected for use in extrapolation is correct, expression (9) below allows us to combine the confidence bound above with the Lipschitz condition, to obtain a conservative test of significance at the  + 1 level using P(C  > f (, Xj ) ≡ g(, A( , Xj ))) > 1 − ,

(9)

and testing ((yj − C )/ ˆ ) > 0 at the 1 level. In the test statistic yj is the observed data at a pre-selected point at condition j, and C is the constant from the Lipschitz condition on g. Sometimes such a constant can be determined from physical theory. Other times C may be estimated, in which case a Bonferroni approach to significance level evaluation results in a test at the  + 1 + 2 ˆ In particular, although the test statistic (yj − C )/ ˆ may be level where 1 − 2 is the level of confidence in an upper bound for C. much greater than 0 at some very high significance level, the true level of the test has to reflect the variation in our estimate of .

1714

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

What this approach does is substitute a computationally simple adjustment of the level of the test for what may be a computationally challenging more accurate estimate of the error variance. As the term “risk orthogonal” implies, this is particularly for use when  and hence C are small. For an example of a case where Cˆ may be estimated, consider the example given in Section 3. The problem considered there is radiation induced darkening of an optical fiber. Ionizing radiation causes the formation of color centers in the glass. In the dose and dose rate range we consider in the example, the induced loss at any given wavelength, measured in decibels (dB) per unit length, is linear in the concentration of these color centers. Since C derives from f (, X) = g(, A( , X)), which is our estimated ˆ , X)) to derive an upper bound for C, ˆ if the example lent itself function of degradation, we could use confidence bounds on g(ˆ , A( to risk orthogonality. Alternatively in the non- risk orthogonal case, the test statistic takes the form:  ⎛     2 ⎝ jf (ˆ ,Xj ) jf (ˆ ,Xi )  ˆ j j

(yj − f (ˆ , Xj )) ⎞.



 −1 jf (ˆ ,Xj ) jf (ˆ ,Xi ) + 1⎠ j j i=1, ...,n i=1, ...,n

(10)

When it can be used, the approach based on risk orthogonality can give quite dramatic results (e.g. LuValle et al., 1998) where the experiments show a signal where none should be detectable. The use of sequential experiments based on (10) has certain interesting problems to overcome. The first term in parenthesis in the denominator of (10) can deflate quite dramatically as the experimental region get closer to the region being extrapolated to, under the null hypothesis of no evanescent processes. On the other hand, further experiments extending toward the region of extrapolation can introduce both a bias into the estimated mean, and inflation to the estimated standard error if evanescent processes exist. In the third section, we discuss conservative estimates of these effects based on a combination of simulation, and approximations using activity maps to estimate the power of prospective experiment sequences in detecting minimally disruptive alternatives. 2.2. Numerical example of a risk orthogonal experiment In the remainder of this section, we extend the example started in the previous section to an example where design of risk orthogonal experiments are nontrivial. The data are generated assuming a uniform distribution of thresholds, a first order kinetic process, an activation energy of 2.25 eV, and a value of  = 1016 Hz. Suppose that we run an experiment at 300 ◦ C for 6 h, with 500 devices, monitoring failure times. When we see that a large number of devices have failed (98%) we run an experiment at 250 ◦ C for 48 h, again with 500 devices. The data is generated using an exponential distribution (a uniform distribution of thresholds across integrated activity for these experiments). The data is plotted on Weibull plots for each simulated experiment in Fig. 5. In the titles, a is the maximum likelihood estimate of the scale parameter, b is the maximum likelihood estimate of the shape parameter, and Pfail is the proportion of devices that failed in the test. Using exponential likelihood calculations based on the observed Weibull shape parameters, the activation energy is estimated as 2.27 eV, with 97.5% two-sided confidence bounds of 2.17--2.35 eV. The 97.5% confidence bounds of log10 () assuming the uniform distribution of thresholds and based on the confidence bounds for the exponential parameter of the first experiment, divided by exp(−Eˆ a /kT) is 15.68--17.39. We can use the joint 95% confidence box generated by the confidence bounds for (Ea , log10 ()) to calculate a 95% upper bound for , the integrated activity under the model, for the experiment run for 1 year at 100 ◦ C. Fig. 6 shows six plots. On the right are the perspective plots, and the left, the image plots for activity maps. The top pair corresponds to Fig. 3, only now it is the activity map showing maximum integrated activity for both the 6 h 300 ◦ C experiment and the 48 h 250 ◦ C experiment vs. life. The middle pair of plots corresponds to Fig. 4, it includes the 1 year 100 ◦ C experiment with the 6 h 300 ◦ C experiment and the 48 h 250 ◦ C experiment in calculating the maximum activity during the experiments then subtracts the maximum activity during life. The bottom plot is the activity map subtracting the integrated activity for the 1 year 100 ◦ C experiment from the maximum integrated activity of the 6 h 300 ◦ C experiment and the 48 h 250 ◦ C experiment. The 95% confidence bound of the signal expected from the 100 ◦ C experiment, , is calculated by calculating the maximum of the integrated activity of the 100 ◦ C experiment over the box on the image plot on the bottom right. The estimated value of  from the activity map is 6 × 10−6 , in this case translating directly to a probability of a single failure in a new device. However, a referee's comments prompted me to examine the plot over only the rectangle with a finer grid including the endpoints, and the value of  is 4 × 10−5 at the extreme corner (Ea = 2.17, log10 () = 17.4). This is shown in Fig. 7. The values on the z-axis are calculated in Fig. 7 are   Ea . z = 1 − exp −t ×  exp − kT Time, t is in seconds,  is in Hertz, Ea is in electron volts, T is temperature on the Kelvin scale, and k is Boltzmann's constant = 8.625 × 10−5 eV/K. Using the same sample size as in the other experiments, the probability of no failures in 500 devices is greater than 0.98 for this value of . Even one failure in this experiment would be a reasonable indication that the simple model of a single Arrhenius

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

300 C Weibull plot a = 1.5 b = 1 Pfail = 0.978

1715

250 C Weibull plot a = 40.57 b = 1.1 Pfail = 0.324

5.000

0.100

0.500

-log (surv)

-log (surv)

0.050

0.050

0.010

0.005 0.005

0.001

0.010

0.100 time

1.000

0.1

0.5 1.0

5.0 10.0

50.0

time

Fig. 5. Weibull plots of data from simulated experiments at 300 and 250 ◦ C.

process with an exponential failure mode is incorrect. However, the total level of the test is 0.07 because of the way  is obtained. To proceed further, we note for example, that with this sample size, there would be strong evidence against the lognormal distribution, because curvature would show in the lognormal plot. However, the conclusion that the shape parameter is 1 may be suspect. To bound epsilon under a Weibull distribution would simply require evaluating the maximum of  over a 1 −  confidence set for the shape parameter, Ea and . Thus for an exponential world, the 100 ◦ C 1 year experiment is conservatively a (0.05, 4 × 10−5 ) risk orthogonal experiment for the observed failure mode for a single Arrhenius process. For a world where lognormal and Weibull distributions are considered to cover the possibilities, the value of either  or  would necessarily increase. 3. Example using activity maps in evaluating power of sequential experiments 3.1. Kinetic model: sequential experiment In this section, we consider a situation where the kinetic model of the degradation mode precludes a risk orthogonal experiment. In particular, this arises in the study of the effect of high energy ionizing radiation on glass containing both pure silica, and boron doped silica components (LuValle et al., 2006). In this case, incremental loss at any given wavelength is proportional to a positive linear combination of the concentrations of different radiation generated defects. Further it is reasonable to assume the defects are generated according to a power law in the dose rate of the radiation (for a monochromatic radiation source, or a source with a single photon energy) and any given class of defect anneals according to a distribution of activation energies in some nth order kinetic process. The simplest form of process involving growth and annealing of radiation effects has the form: d [A1 ] = −(k1new []new )[A1 ] + knew1 [Ai ], dt d[Anew ] = k1new []new [A1 ] − knew1 [Anew ]. dt

(11)

This is the form we assume for the evanescent process. Because any process that involves radiation exposure will cause some darkening, there is no accelerated test that can be made reasonably risk orthogonal, so it is absolutely necessary to calculate the

1716

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

15 1 log10 (nu)

Z

0.5 0 -0.5

5

-1

15

2

5

10

log1

5

0

0 (nu

5

-5

)

0.

1

1.

3

Ea

20

2.

5

0

-10 0.0

-10

0.5

1.0

1.5

2.0

2.5

3.0

2.0

2.5

3.0

2.0

2.5

3.0

Ea 15

Z

0 -0.5

5

-1

15

5

10

log1

5

0 (nu

0 5 )

5

0.

1

1.

2

3

Ea

20

2.

log10 (nu)

1 0.5

5

0

-10 0.0

0.5

1.0

-10

1.5 Ea

15

Z

0 -0.5

5

-1

15

5

10

log1

5

0 (nu

0 )

5

-5

0.

1

1.

2

Ea

20

2.

3

log10 (nu)

1 0.5

5

0

-10 0.0

0.5

1.0

-10

1.5 Ea

Fig. 6. Activity map for risk orthogonal calculation: Top pair of plots, combined 6 h 300 ◦ C experiment and 48 h 250 ◦ C experiment vs. 25 years 25 ◦ C life. Middle pair of plots, combine 1 year 100 ◦ C experiment with other experiments vs. 25 years 25 ◦ C life. Bottom pair of plots, combined 6 h 300 ◦ C experiment and 48 h 250 ◦ C experiment vs. 1 year 100 ◦ C experiment, with overlaid confidence box for(log10 (), Ea ).

variance for prediction from the chosen base model. In this particular case the base model has the following form: ⎞ ⎛ 11 21 11 21     d[A1 ] 1 2 ⎠ ⎝ k1i [] + k1i [] ki1 [Ai ]n1 + n−1 ki1 [Ai ]n2 , [A1 ] + n−1 = − 1 2 dt i=2

i=12

d[Ai ] = k1i nj []j [A1 ] − ki1 [Ai ]nj , dt

i=2

i = 2, . . . , 21, j = 1, 2,

i=12

(12)

where two populations of defects were identified, each with a distribution estimated by a mixture of beta distributions, and n1 = n2 = 2. This is slightly different from the analysis of this data previously published (LuValle et al., 2006). The reasons for the form chosen in (12) are the following: (1) Experience with radiation induced defects in glass has shown that each population of defects anneals according to a distribution of activation energies, but has a single value of . Physical arguments can be made for this based on explicit form of  as a difference in entropies between reactants and the activated state.

Integratedactivity ∗1e5

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

1717

5 4 3 2 1 0 17 16 .5 10 (nu )

2.3

log

16

2.25 2.2

Ea

Fig. 7. Integrated activity for 1 year 100 ◦ C experiment over the confidence box for(Ea , log10 ()).

(2) Graphical analysis of this data (LuValle et al., 2006) at multiple wavelengths of light over time indicate at least two populations of defects. (3) This model had the lowest sum of squared residuals of the > 30 models representing a sequential model, and various combinations of kinetic order of the annealing (LuValle et al., 2006). The selection of second-order kinetics rather than first-order kinetics for both anneals comes from using higher resolution in time than was used in the previous analysis. In this model [] is dose rate, the k1i , i = 2, . . . , 11 is distributed evenly between 0 and 1 to approximate a mixture of beta densities: m(x; p, 11 , 12 , 13 , 14 ) = p  kij = j exp −

Eaij kT

 ,

j=

(11 + 12 ) 11 −1 (13 + 14 ) 13 −1 (1 − x)12 −1 + (1 − p) (1 − x)14 −1 x x (11 ) (12 ) (13 ) (14 )

1 2

(13)

ij = 2, . . . , 11 ij = 12, . . . , 21

while Eaij = Eajlow + (Eajhigh − Eajlow )

ij − min(ij − 1) . 11

(14)

Expression (14) shows the way the rate constants for annealing change with i and j. In particular, the Arrhenius premultiplier for the anneal,  is assumed constant for each population, and the two separate populations are indexed by j. The index i (denoted ij above) is essentially nested in j, with each value of i representing a subpopulation of defects with an identical activation energy. The model is a nonlinear system of differential equations (in contrast with compartment models for example), so the fitting (and covariance calculation using numerical derivatives) takes a long time. Thus, any scheme to reduce computational time is important. For our evanescent process, (11), we need to specify specific values of k1new , new , Eanew , new where knew1 = new exp(−Eanew /kT). Once the evanescent process is specified, the second task is finding the values of these parameters that will result in significant activity over operating conditions but almost no activity over the course of the experiments to date. For the calculations below to do this, the integrated activity during the (actual and proposed) experiments was calculated over the four dimensional array of parameters given by log10 (k1new ) = (−10, −9, . . . , 20),

new = (0.2, 0.4, . . . , 3.8, 4), log10 (new ) = (−9.99, −6.66, . . . , 16.65, 19.98), Eanew = (0.3, 0.6, . . . , 2.7, 3.0). The integrated activity at each point gives us some very useful information. Loss in an optical fiber takes a particularly simple form. The loss in dB is proportional to the concentration of loss causing defects in the fiber, which is proportional to

1718

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

the integrated activity if we assume that the defects do not interact in some nonlinear manner to produce higher loss as the concentration increases. 3.2. Theoretical and computational approximations To see how this helps, assume that there is a value of loss, X0crit , such that if the loss rises above that value, the device is considered failed, so predicting increase to such a loss means the device will not be sent to market (assume it is not field replaceable). Then by the linearity of loss with normalized integrated activity, at an experimental condition conditioned on a set of parameter values we can define Tcrit.i = (X0crit − Xpred0 )([Anew ]i. max /[Anew ]life. max ), where [Anew ]∗. max is simply proportional to the maximum integrated extent during experiment i (life), and Xpred0 is the predicted loss under the base model. Tcrit.i is the critical increased loss we would have to be able to detect at experiment i under the assumed parameter values in order to identify the evanescent process. Using delta method asymptotics, and assuming Gaussian measurement error, an estimate of the probability of detecting this critical value using a simple one-sided statistical test of level  is  1 − z −

Tcrit.i spred,i−1,i

.

(15)



ˆ 2pred.i−1,i + ( ˆ 2error /nexp t,i ), ˆ 2pred.i−1,i is the standard error of the mean prediction at the ith experiment using 2 all experiments up to and including the i-1st, ˆ error is the standard deviation of the residuals, nexp t,i is the sample size of the ith experiment and is the standard Gaussian CDF. This estimate is clearly correct for testing a single point of the curve if the

Here spred.i−1.i =

process is truly evanescent, for all the experiments being used to estimate the standard errors and Xpred0 . However if there are a sequence of experiments in which the processes become less evanescent, but not large enough to be statistically significant, they will tend to inflate the standard error terms and deflate Tcrit. exp t . To estimate the effect of deflation, we calculate numerically the derivative of the prediction with respect to a small change in the expected value of the response at a new experiment dXpred0 /de and re-estimate ⎛

⎛

⎛

T˜ crit.i = ⎝X0crit − ⎝Xpred0 + ⎝

 j
⎞⎞⎞  dXpred0 ⎠⎠⎠ × [Anew ]i. max . Tcrit.j dj [Anew ]life. max

(16)

To bound the effect of inflation on the standard error for each added experiment in the sequence the following adjustment is conjectured to be conservative:

ˆ 2error.adj.i

   2 2 2 , (z s (nbase ) ∗ ˆ error + j




 2  ˆ error.adj.i−1

ˆ 2error 2  ˆ

pred,i−1,i + spred.adj.i−1,i = , nexp t

ˆ 2error

(17)

where nbase is the number of experiments used in identifying the original model. The term  j
  2 min Tcrit.j , (z spred.adj,j−1,j )2 z −

Tcrit.j spred.adj,j−1,j

,

2 is an upper bound on the squared bias contributed by the evanescent process. The term Tcrit,j is the squared bias assuming no testing at experiment j and we just add in the results of the new experiment uncritically to improve our variance estimate for further experiments. The term (z spred.adj,j−1,j )2 (z − (Tcrit.j /spred.adj,j−1,j )) is the square of the largest value that our t test will

have and still not reject the alternative model times the probability that the observed test statistic is smaller than or equal to that 2 value. The sums in the definition of ˆ error.adj.i and in Eq. (16) are over the experiments preceding experiment i used to predict the value observed during experiment i. The sum in the denominator of the first line of (17) counts the number of experiments included in evaluating the variance. The two conflicting effects, the reduction of the prediction variance under the base model by experiments approaching operating conditions, and the potential deflation of the signal and inflation of error variance that can arise when an evanescent process exists, makes for an interesting problem in design of inference protocols. Ideally one would use a pattern of pooling that would optimize the signal to noise ratio for testing. By using the activity maps as above to explore the power of a given sequences of tests, it is possible to pre-choose a pooling/testing strategy adapted roughly to the given problem. With significant computing resources, one could select a large number, N, of proposed experiments, and explore the power of all 2N − 1 patterns of

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

1719

Table 1 New experiments to explore Dose rate (rad/s)

Time (h)

T (◦ C)

Finished experiments A B C

2.2 110 110

12.6 0.25 2.5

25 25 25

New experiment # 1 2 3

2.2 2 × 10−2 2 × 10−2

31.5 720 5000

5 25 25

Table 2 Power calculations for proposed new experiment sequence Base experiments

1, 0.01 level

1, 0.03 level

2, 0.01 level

2, 0.03 level

3, 0.01 level

3, 0.03 level

ABC ABC1 ABC12

0.02 X X

0.05 X X

0.08 0.19 (0.20) X

0.13 0.25 (0.27) X

0.21 0.45 (0.46) 0.38 (0.55)

0.27 0.55 (0.57) 0.48 (0.66)

experiments, with associated multiple patterns of pooling and testing. Interestingly, some initial simulations using a linear model analog to evanescent processes indicated that nearly optimal power could be obtained by testing a pooled set of experiments vs one new experiment at a time. A simplified approach that takes advantage of this observation is the following. To compute approximate power: (i) for the n original experiments used to explore the data, and m new experiments; (ii) use the base model from the original experiments, use (15) to construct estimates of the power of a simple t test at the point of maximum signal for each of the m new experiments using the Tcrit values calculated from the activity maps; (iii) add the experiment point closest to the n original experiments to the pool estimating the base model (using data simulated from the base model). Then use (15)--(17) to construct estimates of the power of a simple t test at the point of maximum signal for each of the (m − 1) remaining experiments, calculating the values of the various Tcrit type terms using the activity maps; (iv) repeat (iii) until only one of the experiments is tested, and the remaining m − 1 are pooled. The variance estimates (assuming the base model) require large amounts of computing time, so I have found it reasonable to separate the variance and dXpred0 /de terms (both of which can be calculated using the base model) for each pooling as a separate run. Then the results of those calculations can be combined in a list, and used against a list from the activity map calculations to arrive at the final power calculations. An important point to recall in the power calculation using any Tcrit type terms, is that the inference sequence is followed separately for each point in the parameter space in the part of the activity map where evanescent processes are indicated. The calculation of a single number to represent power thus requires some accepted scheme of combining parameter values over the section of space where the activity map was negative for the base set of experiments. The approach taken here is to assume a uniform prior over this region, although, use of any easily specified prior would be a simple extension. 3.3. Numerical results Table 1 shows the base sequence of experiments (actually performed) and a proposed new sequence. As in LuValle et al. (2006), each experiment also includes a 500 h annealing period at room temperature, and we are extrapolating to a 10-year, 100-krad dose at 25 ◦ C. The predicted end of life loss under those conditions is ∼0.2 dB/km. The minimal loss causing a system problem is assumed here to be 1 dB/km at the wavelength of interest. The approximate power calculations are based on applying the scheme for different test levels for individual tests and are shown in Table 2. The numbers outside of parentheses are the integrated powers for the individual tests. The numbers in parentheses represent the integrated power for the whole procedure including all tests to the left. So if we plan to do three tests at the 0.03 level: 1. test the results of experiment 1 vs. the extrapolation to the conditions of experiment 1 using the original data to estimate the base model and prediction error; 2. test the results of experiment 2 vs. the extrapolation to the conditions of experiment 2 using the original data plus the data from experiment 1 to estimate the base model and prediction error;

1720

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

3. test the results of experiment 3 vs. the extrapolation to the conditions of experiment 3 using the original data plus the data from experiment 1 plus data from experiment 2 to estimate the base model and prediction error. The power of the joint 0.09 test, calculated in sequence at each point in the volume of uncertainty using formula 15--17, and then integrated over the volume, is 0.66. Note that as we proceed down the last column, looking at the number outside parentheses, the integrated power is maximized pooling only experiment 1. Adding experiment 2 to the estimate distorts the test statistic enough under the alternative that the integrated power of that individual test is now reduced. This reduction in calculated power is because the addition of new data, using the sequential testing and estimation scheme outlined above, is adding bias to the test statistic that reduces sensitivity. An interesting conjecture is that this pattern of diminishing returns on power for detecting departures from the model used for extrapolation is fundamental to the extrapolation problem with an unknown model, and not just a result of the simple testing formulation used here for computational exploration. Although these calculations are based on several approximations, they agree qualitatively with the results of simulations made using a linear model analog to evanescent processes, where the evanescent process was emulated by a change in slope at some point between the region where data was gathered and the region being extrapolated to, and change in slope was calculated to give a particular increase at the point being extrapolated to. There are several reasons one might expect that these power calculations are conservative: (A) an actual test would be along the whole observed loss growth and annealing path; (B) the calculations in (15)--(17) are designed to be conservative; (C) the alternative is the minimal alternative which would cause a problem in the system; any worse disruptive alternative would be easier to detect. For example, if the loss at operating conditions is actually 2 dB (total loss) then the number in parenthesis on the lower right becomes 0.89. It is, however possible that the use of linearization, the dXpred0 /de terms, to calculate bias will sometimes result in the power estimate being overly optimistic.

4. Summary and open questions It is plain from the Stone--Weierstrass theorem that any extrapolation based on purely statistical methods of assessing fit, and in particular using no input from subject matter science is bound to have a serious problems with confounding. What is interesting is that even when rigorous attention is paid to using all subject matter information available, this confounding problem may still exist. Further, such confounding has evidenced itself in real life with problems that have had serious impact on the telecommunications and electronics industries, as described in the introduction. Here and in LuValle (2007) a new tool, activity maps, have been introduced which allow precise identification of potential inference disrupting, confounded alternative models. In LuValle (2007) activity maps are defined and explored for a number of simple cases of kinetic models, and their relationship to demarcation maps is elucidated. In this paper, they are integrated into statistical theory to provide a tool for approximately evaluating the power in designing experiments to explore for potential disruptive alternatives. Two uses of activity maps, one for precise identification of risk orthogonal experiments, and the potential signal they provide, and the other for assessment of the power of a sequence of experiments approaching the region being extrapolated too, have been discussed, and demonstrated with examples. The first with a simulation (although a rough version based on demarcation maps, where the experiment was performed and identified an evanescent process has been described; LuValle et al., 1998), the second with a real example. In the second, where risk orthogonal experiments are not possible, an interesting trade off in sequential design for extrapolation was found. In particular, pooling data into the base model as the independent variables governing the phenomena of interest approach the condition being extrapolated to can result in either of two conflicting phenomena: (1) The prediction variance for subsequent data even closer to the conditions of interest can be substantially reduced if no evanescent process is active in the data being pooled. (2) The error variance can be inflated, and the prediction of the base model substantially biased if an evanescent process is active, but not active enough to be detected statistically. A simple approach to balancing these conflicting phenomena, using a sequence of poolings and tests, reminiscent of setting a baseline for a quality control chart, then extending that base set, has been proposed, and implemented in the real example. However, more sophisticated approaches will almost surely improve the statistical power of the approach. For example, perhaps the theory of Pascual (2006) for making designs robust to structural model misspecification, based on White's (1982) quasi maximum likelihood estimation could be applied here. The setup is similar to the problem of watching dynamic systems as the various variables controlling the systems change, (e.g., hurricane generation as ocean temperatures warm, or ecosystem change as carbon dioxide levels rise) so perhaps similar calculations can be used to help plan observational studies designed to test particular models in such systems.

M.J. LuValle / Journal of Statistical Planning and Inference 139 (2009) 1707 -- 1721

1721

Finally, one particular point needs to be noted here. This design theory is quite different from most of that applied to accelerated testing. In particular, we assume that there is economic reason to protect ourselves against problems, sufficient that we are willing to extend the design of accelerated tests to quite long time frames (even if only to give pre-warning of potential problems in the field). This assumption is (barely) correct in some parts of the telecommunications industry. In health, defense, and aerospace industries where loss has higher cost, it should be more easily justified. References Augis, J., Denure, D., LuValle, M., Mitchell, J., Pinnel, M., Welsher, T., 1989. A humidity threshold for conductive anodic filaments in epoxy glass printed wiring boards. In: Proceedings of the Third International SAMPE Electronics Conference, pp. 1023--1030. Bagdonavicius, V., Nikulin, M., 2002. Accelerated Life Models: Modeling and Statistical Analysis. CRC Press, Boca Raton. Chermoff, H., 1972. Sequential Analysis and Optimal Design. SIAM, Philadelphia. Elfving, G., 1952. Optimal allocation in linear regression theory. Ann. Math. Statist. 23, 255--262. Glasstone, S., Laidler, K., Erying, H., 1941. The Theory of Rate Processes. McGraw-Hill, New York. Klinger, D., 1992. Failure time and rate constant of degradation, an argument for the inverse relation. Microelectron. Reliability 32, 987--994. Laidler, K.J., 1987. Chemical Kinetics. third ed. Harper & Row, New York, NY. Lenat, D., Brown, J., 1984. Why AM and EURISKO appear to work. Artif. Intelligence 23 (3), 269--294. LuValle, M., 1999. An approximate kinetic theory for accelerated testing. IIE Trans. Quality Reliability 31 (12), 1147--1156. LuValle, M., 2004a. Some Bayesian experimental design theory for risk reduction in extrapolation. J. Risk Anal. V24 (5), 1249--1259. LuValle, M., 2004b. A mathematical framework for identifying neighborhoods of a certain class of nonlinear dynamic models for extrapolating observed degradation phenomena. In: Antonov, V., Huber, C., Nikulin, M., Polischook, V. (Eds.), Proceedings of the Longevity Aging and Degradation Models in Reliability, Public Health, Medicine and Biology (LAD 2004), St. Petersburg State Politechnical University, St. Petersburg, pp. 155--166. LuValle, M., 2007. Identifying mechanisms that highly accelerated tests miss. IEEE Trans. Reliability 56 (2), 349--359. LuValle, M., Hines, L., 1992. Using step stress to explore the kinetics of failure. Quality Reliability Eng. Internat. 8, 361--369. LuValle, M., Welsher, T., Svoboda, K., 1988a. Acceleration transforms and statistical kinetic models. J. Statist. Phys. 52, 311--320. LuValle, M., Copeland, L., Kannan, S., Judkins, J., Lemaire, P., 1998b. A strategy for extrapolation in accelerated testing. Bell Labs Tech. J. 3 (3), 139--147. LuValle, M., Mrotek, J., Copeland, L., Lefevre, B., Brown, G., Throm, R., 2002. Integrating computational models with experimental data in reliability. In: Matthewson, M., Kurkjian, C. (Eds.), Proceedings of the SPIE, vol. 4215, pp. 52--63. LuValle, M., Lefevre, B., Kannan, S., 2004. The Design and Analysis of Accelerated Tests for Mission Critical Reliability. CRC Press, Boca Raton. LuValle, M., Friebele, E., Dimarcello, F., Monberg, E., Wasserman, L., Wisk, P., Yan, M., Birch, E., 2006. Radiation induced loss predictions for pure silica core polarization maintaining fiber, Photonics Europe. In: Limberger, H., Matthewson, M. (Eds.), Proceedings of the SPIE, vol. 6193, 61930J. Matthewson, M.J., 1998. Chemical kinetics models for the fatigue behavior of fused optical silica fiber. Mater. Res. Soc. Symp. Proc. 531, 143--153. Meeker, W.Q., 1987. Limited failure population life tests: application to integrated circuit reliability. Technometrics 29, 51--65. Meeker, W.Q., Escobar, L., 1998. Statistical Methods for Reliability Data. Wiley, New York. Meeker, W., Hahn, G., 1985. How to plan and accelerated life test-some practical guidelines. In: ASQC Basic References in Quality Control: Statistical Techniques, vol. 10. Meeker, W.Q., LuValle, M.J., 1995. An accelerated life test model based on reliability kinetics. Technometrics, 133--148. Meeker, W.Q., Nelson, W., 1975. Optimum accelerated, life tests for Weibull and extreme value distributions and censored data. IEEE Trans. Reliability R-24, 321. Nelson, W., 1990. Accelerated Testing: Statistical Models, Test Plans, and Data Analysis. Wiley, New York. Pascual, F.G., 2006. Accelerated life tests robust to misspecification of the stress--life relationship. Technometrics 48 (1), 11--25. Primak, W., 1955. The kinetics of processes distributed in activation energy. Phys. Rev. 100 (6), 1677--1689. Turner, T., Wendel, K., 1985. The influence of stress on aluminum conductor life. In: Proceedings of the International Reliability Physics Symposium, pp. 142--147. Viertl, R., 1988. Statistical Methods in Accelerated Life Testing. Vandenhoeck, Ruperecht, Gottingen. Welsher, T., Mitchell, J., Lando, D., 1980. CAF in composite printed-circuit substrates: characterization, modeling, and a resistent material. In: Proceedings of the International Reliability Physics Symposium, pp. 235--237. White, H., 1982. Maximum likelihood estimation of misspecified models. Econometrica 50, 1--25. Yost, F.G., Amos, D.E., Romig, A.D., 1989. Stress driven diffusive voiding of aluminum conductor lines. In: Proceedings of the International Reliability Physics Symposium, pp. 193--201.